{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 2: Introduction to Predictive Modeling\n",
    "\n",
    "> Ignorance is preferable to error and he is less remote from the truth who believes nothing than he who believes what is wrong.\n",
    "Thomas Jefferson (1781)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objectives\n",
    "+ To define predictive modeling.\n",
    "+ To introduce the idea of structural causal models and their graphical representation.\n",
    "+ To tell the difference between aleatory and epistemic uncertainties.\n",
    "+ To introduce the uncertainty propagation problem.\n",
    "+ To introduce the model calibration problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Readings\n",
    "\n",
    "+ These notes.\n",
    "\n",
    "+ [Oden, Moser, Ghattas, Computer Predictions with Quantified Uncertainty, Part I](https://archive.siam.org/pdf/news/1842.pdf)\n",
    "\n",
    "+ [Oden, Moser, Ghattas, Computer Predictions with Quantified Uncertainty, Part II](https://archive.siam.org/pdf/news/1857.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the things we need to plot\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "sns.set_context('talk')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predictive Modeling\n",
    "\n",
    "> Reports that say that something hasn't happened are always interesting to me, because as we know, there are known knowns; there are things we know we know. We also know there are known unknowns; that is to say we know there are some things we do not know. But there are also unknown unknowns—the ones we don't know we don't know. And if one looks throughout the history of our country and other free countries, it is the latter category that tend to be the difficult ones. Donald Rumsfeld, United States Secretary of Defense, [DoD news briefing, February 12, 2002](https://archive.defense.gov/Transcripts/Transcript.aspx?TranscriptID=2636).\n",
    "\n",
    "*Predictive modeling* is the process of describing our state of knowledge about known unknowns in order to make informed decisions.\n",
    "This is the scope of this class.\n",
    "\n",
    "Unfortunately, there is no automated way for turning unknown unknowns to known unknowns.\n",
    "This is currently done manually as it requires common sense and human intuition.\n",
    "Automating this process seems to require the ability to perform induction on open-ended problems and it may require general artificial intelligence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Structural causal models\n",
    "\n",
    "A *causal model* is a model that attempts to capture the mechanisms that govern a given phenomenon.\n",
    "We will use the language of *structural causal models* (SCM), developed by the computer scientist Judea Pearl, to formalize the concept.\n",
    "A structural causal model is a collection of three things:\n",
    "+ A set of variables. These are variables that our model is trying to explain (endogenous), but also other variables that may just be needed (exogenous).\n",
    "+ A set of functions that give values to each variable based on the values of all other variables.\n",
    "\n",
    "Most physical and engineering models are causal models.\n",
    "\n",
    "### Example: Asthma model (J. Pearl)\n",
    "\n",
    "Suppose that we are trying to study the causal relationships between a treatment $X$ and lung function $Y$ for individuals who suffer from asthma.\n",
    "However, it is plausible that $Y$ also depends on the air pollution levels $Z$.\n",
    "The final ingredient is the set of function that connects $X$ and $Z$ to $Y$.\n",
    "$$\n",
    "Y = f(X, Z).\n",
    "$$\n",
    "\n",
    "### Graphical representation of causal models\n",
    "Every SCM is corresponds to a *graphical causal model*.\n",
    "These are usually *directed acyclic graphs* (DAGs).\n",
    "These can be read trivially from the SCM form.\n",
    "Let's look at an example.\n",
    "\n",
    "### Example: Asthma model - Graphical causal model\n",
    "Here I am representing each variable with a node.\n",
    "The node at the beginning on an arrow is the direct cause of the node at the end of the arrow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: Asthma Pages: 1 -->\n",
       "<svg width=\"276pt\" height=\"116pt\"\n",
       " viewBox=\"0.00 0.00 276.07 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
       "<title>Asthma</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-112 272.07,-112 272.07,4 -4,4\"/>\n",
       "<!-- X -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>X</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"57.99\" cy=\"-90\" rx=\"57.97\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"57.99\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">X (treatment)</text>\n",
       "</g>\n",
       "<!-- Y -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>Y</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"128.99\" cy=\"-18\" rx=\"49.14\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"128.99\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\">Y (asthma)</text>\n",
       "</g>\n",
       "<!-- X&#45;&gt;Y -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>X&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M74.81,-72.41C83.83,-63.52 95.1,-52.41 105.02,-42.63\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"107.62,-44.98 112.28,-35.47 102.7,-40 107.62,-44.98\"/>\n",
       "</g>\n",
       "<!-- Z -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>Z</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"200.99\" cy=\"-90\" rx=\"67.17\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"200.99\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">Z (air pollution)</text>\n",
       "</g>\n",
       "<!-- Z&#45;&gt;Y -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>Z&#45;&gt;Y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M183.92,-72.41C174.64,-63.39 163.02,-52.09 152.86,-42.21\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"155.04,-39.45 145.43,-34.99 150.16,-44.47 155.04,-39.45\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a25b59ed0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from graphviz import Digraph # TODO add pygraphviz to dependencies\n",
    "g = Digraph('Asthma')\n",
    "g.node('X', label='X (treatment)')\n",
    "g.node('Y', label='Y (asthma)')\n",
    "g.node('Z', label='Z (air pollution)')\n",
    "g.edge('X','Y')\n",
    "g.edge('Z', 'Y')\n",
    "#g.render('asthma_graph', format='png') # Uncomment the line if you want to save the figure\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Types of Uncertainty\n",
    "In general, we are uncertain about something if we don't know much about it.\n",
    "In particular, we can be uncertain about:\n",
    "+ the value of a model parameter;\n",
    "+ the initial conditions of an ordinary differential equations;\n",
    "+ the boundary conditions of a partial differential equation;\n",
    "+ the value of an experimental measurment we are about to perform;\n",
    "+ the mathematical form of a model;\n",
    "+ etc.\n",
    "\n",
    "Uncertainty may be *aleatory* or *epistemic*.\n",
    "\n",
    "+ Aleatory uncertainty is associated with inherent system randomness. \n",
    "+ Epistemic uncertainty is associated with lack of knowledge.\n",
    "\n",
    "There is a long philosophical debate about this distinction.\n",
    "We are going to ignore it.\n",
    "The instructors view is that common sense and **probability theory are sufficient to describe both uncertainties**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Driving a trailer on a rough road (1/3)\n",
    "\n",
    "The following example is modified from Purdue's Basic Mechanics II Lecture Book.\n",
    "In the figure you see a trailer of mass $m$ moving on a rough road with speed $v$.\n",
    "The suspension spring constant is $k$.\n",
    "We are intersted in the vibration amplitude $X$.\n",
    "\n",
    "![Trailer](trailer.png)\n",
    "\n",
    "We do not know (yet) how we can model a true road, so let us assume that the road surface is sinusoidal with amplitude $y_0$ and \"wavelength\" $L$.\n",
    "Doing a [little bit dynamics](https://www.youtube.com/watch?v=d4OKBqr_aYQ&feature=youtu.be), shows that the amplitude of the suspension oscilation is:\n",
    "$$\n",
    "X = \\left|\\frac{ky_0}{k-m\\omega^2}\\right|,\n",
    "$$\n",
    "where the angular velocity is:\n",
    "$$\n",
    "\\omega = \\frac{2\\pi v}{L}.\n",
    "$$\n",
    "Let's draw the causal graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: Trailer Pages: 1 -->\n",
       "<svg width=\"314pt\" height=\"188pt\"\n",
       " viewBox=\"0.00 0.00 314.00 188.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 184)\">\n",
       "<title>Trailer</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-184 310,-184 310,4 -4,4\"/>\n",
       "<!-- k -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>k</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">k</text>\n",
       "</g>\n",
       "<!-- X -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>X</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"135\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"135\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\">X</text>\n",
       "</g>\n",
       "<!-- k&#45;&gt;X -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>k&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M45.81,-76.81C63,-65.67 88.62,-49.06 107.99,-36.5\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"109.92,-39.43 116.4,-31.05 106.11,-33.56 109.92,-39.43\"/>\n",
       "</g>\n",
       "<!-- m -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>m</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">m</text>\n",
       "</g>\n",
       "<!-- m&#45;&gt;X -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>m&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M107.35,-72.76C111.71,-64.28 117.15,-53.71 122.04,-44.2\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"125.23,-45.64 126.7,-35.15 119.01,-42.44 125.23,-45.64\"/>\n",
       "</g>\n",
       "<!-- y0 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>y0</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"164\" y=\"-86.8\" font-family=\"Times,serif\" font-size=\"14.00\">y</text>\n",
       "<text text-anchor=\"start\" x=\"171\" y=\"-86.8\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\">0</text>\n",
       "</g>\n",
       "<!-- y0&#45;&gt;X -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>y0&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M162.65,-72.76C158.29,-64.28 152.85,-53.71 147.96,-44.2\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"150.99,-42.44 143.3,-35.15 144.77,-45.64 150.99,-42.44\"/>\n",
       "</g>\n",
       "<!-- omega -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>omega</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"243\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"237.94\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">ω</text>\n",
       "</g>\n",
       "<!-- omega&#45;&gt;X -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>omega&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M224.19,-76.81C207,-65.67 181.38,-49.06 162.01,-36.5\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"163.89,-33.56 153.6,-31.05 160.08,-39.43 163.89,-33.56\"/>\n",
       "</g>\n",
       "<!-- v -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>v</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"207\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"207\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\">v</text>\n",
       "</g>\n",
       "<!-- v&#45;&gt;omega -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>v&#45;&gt;omega</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M215.35,-144.76C219.71,-136.28 225.15,-125.71 230.04,-116.2\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"233.23,-117.64 234.7,-107.15 227.01,-114.44 233.23,-117.64\"/>\n",
       "</g>\n",
       "<!-- L -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>L</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"279\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"279\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\">L</text>\n",
       "</g>\n",
       "<!-- L&#45;&gt;omega -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>L&#45;&gt;omega</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M270.65,-144.76C266.29,-136.28 260.85,-125.71 255.96,-116.2\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"258.99,-114.44 251.3,-107.15 252.77,-117.64 258.99,-114.44\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a25b8b750>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = Digraph('Trailer')\n",
    "g.node('k')\n",
    "g.node('m')\n",
    "g.node('y0', label='<y<sub>0</sub>>')\n",
    "g.node('omega', label='<&omega;>')\n",
    "g.node('v')\n",
    "g.node('L')\n",
    "g.node('X')\n",
    "g.edge('v', 'omega')\n",
    "g.edge('L', 'omega')\n",
    "g.edge('y0', 'X')\n",
    "g.edge('omega', 'X')\n",
    "g.edge('k', 'X')\n",
    "g.edge('m', 'X')\n",
    "g.render('trailer_g', format='png')\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "+ Which variables could be \"known unknowns?\"\n",
    "+ Which of these are aleatory and which epistemic?\n",
    "+ How can you reduce the epistemic uncertainty of some of these variables?\n",
    "+ What are some \"unkown unknowns\" that you could turn into \"known unknowns?\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The uncertainty propagation problem\n",
    "\n",
    "As we mentioned earlier, both aleatory and epistemic uncertainties can be described using probability theory.\n",
    "So, one of the first steps of predictive modeling is to come up with probability statementes for all uncertain variables.\n",
    "However, this is also one of the most difficult problems...\n",
    "So, let's assume that some has already done it for us.\n",
    "The next step is to propagate this uncertainty through the causal model to characterize our uncertainty about a quantity of interest.\n",
    "Let us do both using the trailer example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving uncertainty propagation problems\n",
    "\n",
    "The simplest way to solve the uncertainty propagation problem is via sampling.\n",
    "This is known as the *Monte Carlo* method. \n",
    "It was invented in Los Alamos during the Manhatan project.\n",
    "We will study the Monte Carlo method extensively.\n",
    "For now, let's look at a simple example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Driving a trailer on a rough road (2/3)\n",
    "\n",
    "To make this more precise, assume that we are the manufacturer of the trailer.\n",
    "Let's quantify our state of knowledge about all the parameters of this model using a little bit of common sense.\n",
    "\n",
    "| Variable | Type | Values| \n",
    "|:---------|:--------------|:--------|\n",
    "| $k$ | Manufacturing uncertainty | [159,999, 160,001] N/m |\n",
    "| $v$ | Operating condition | [80, 150] km/hour |\n",
    "| $m$ | Loading condition | [100, 200] kg|\n",
    "| $y_0$ | Road condition | [0, 100] mm|\n",
    "| $L$ | Road condition | [1, 2] m |\n",
    "\n",
    "Not being able to come up with more precise information (or any data) we would consider any value within this intervals as equally likely.\n",
    "Now, let's write some code to see how this uncertainty affects the angular velocity $\\omega$ and the amplitude $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for numerical arrays and linear algebra:\n",
    "import numpy as np \n",
    "# The number of samples we wish to take\n",
    "num_samples = 1000\n",
    "# Two arrays in which to store the samples we take\n",
    "Xs = np.ndarray((num_samples, )) # To store the samples\n",
    "omegas = np.ndarray((num_samples, ))\n",
    "for i in range(num_samples):\n",
    "    k = 160000.0 + np.random.rand() # np.random.rand() samples a number uniformly between 0 and 1\n",
    "    m = 100.0 + (200.0 - 100.0) * np.random.rand() # Here we sample a random number in [100, 200]\n",
    "    y0 = 100 * np.random.rand() * 1e-3 # Turning it to m\n",
    "    v = (80.0 + (150.0 - 80.0) * np.random.rand()) * 1e3 / 3600.0 # Turning it to m/s\n",
    "    lam = 1.0 + (2.0 - 1.0) * np.random.rand()\n",
    "    omega = 2.0 * np.pi * v / lam\n",
    "    X = np.abs(k * y0 / (k - m * omega ** 2))\n",
    "    omegas[i] = omega\n",
    "    Xs[i] = X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbIAAAFdCAYAAABmYehLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3debgcVbn2/+8NCIQgCYkEVFQUDGF6FUUFX2ZB5imIHgQSECcUOUfhByoRAV9QEEQ5iAcBSeSoOEBANChhCHIUFYQjQyRAAAHDGMhmCmF6fn+s1UmlUz3V7s7em9yf6+qr965atWpVZaefXqvWoIjAzMxsqFpuoAtgZmbWHw5kZmY2pDmQmZnZkOZAZmZmQ5oDmZmZDWkOZGZmNqQ5kJmZ2ZDmQGZmZkPaClUOkjQW2AgYAwTwOHB7RNzdxbKZmZm1pHZn9pC0AfBZYD9gzdrm/F7L5FHgF8A5EfGPLpbTzMysVMtAJmld4BRgH2A+cD1wAzAbmEsKZqOA9YDNga2AYcAlwDERcW+vCm9mZtZOIFsA3AZ8D7gkIp5rkX448BHgCGCjiFi5S2U1MzNbQjuBbK+IuKxS5v041szMrB1tPyMzMzMbjNz93szMhjQHMjMzG9K6GsgkvUXSNZKu7ma+ZmZmjVQaEN3EKsC2LBpXZmZm1lPdDmSzgbd3OU8zM7OG3GvRzMyGNHf2MDOzIa1l06KkCVUyjogfVznOzMysE+3M7PEqqfOGmiZcXETE8v0pmJmZWTva6eyxXc9LYWZmVpE7e5iZ2ZDWr84eklaS9GZJK3arQGZmZp2oFMgkvUfSNcAzwAPAlnn7GElXS9qhi2U0MzNrqONAJundpMU11wUW65kYEY+RFtWc2JXSmZmZtVClRnYiMAfYCPgyS/ZmvBp4fz/LZWZm1pYqgWwr4NyIeJbyORUfAN7Ur1KZmZm1qUogWxnoa7J/tYplMTMz61iVQDYbeG+T/dsDM6sVx8zMrDNVAtlPgYPqeiYGgKQjgZ2BC7tQNjMzs5Y6HhCdx4z9HtgauBMYB9wGrAGsBUwHdo2IV7tbVDMzsyV1XCOLiBeBHYGjgPnAC8BY4AngaGB3BzEzM1taqtTIlo+IV3pUHjMzs45UeUY2R9J3JL2r66UxMzPrUJUa2Q3AB0gdPG4HJgM/jYhHu146MzOzFirNfi9pPdI0VAcA6wAvA1cCU4BfR8SCLpbRzMysoX4v4yJpG2ACsC9pMPQ84BcR8dn+F8/MzKy5rq1HJmllUg3tdOD1XiHazMyWhnZWiG5J0vakWtl4YFVgbjfyNTMza6VyIJM0jhS8DgDWJj0nu4L0nOw3XSmdmZlZC1V6LR5OCmDvJS3hcgspeP00Ip7oegnNzMyaqBLIXgUeBn4CTImIO3pRMDMzs3ZUaVrcFbjS01CZmdlg0LVei2ZmZgOhyhRVZmZmg4YDmZmZDWkOZGZmNqQ5kJmZ2ZDmQGZmZkOaA5mZmQ1pDmRmZjakdTWQSXqLpGskXd3NfM3MzBrpyuz3BasA25JWjzYzM+u5bgey2cDbu5ynmZlZQ56iyszMhjR39jAzsyGt46ZFSde0SBLAfOAB4ErgsnC1z8zMeqTKemT3A8OANfKmefl9ZH5/nFTTG00Kan8EdomI5/pbWDMzs3pVmha3BZ4Hvg2sGRGjImIUsCZwGvAcsBnwBuB0YEvguK6U1szMrE6VGtlU4LmIOLDB/v8GVomI8fn3y4ENImK9/hbWzMysXpUa2XbA9U32/w+wfeH3q4C1K5zHzMyspSrjyASMa7J/XE5T8wqp84cBkl4mfYF4eqDLYmY2hKwGvBoRS8StKoHsKuAwSX+JiIuKOyTtD3wW+E1h82bA/RXO81q1HKARI0aMGOiCmJkNFX19fdCgFbHKM7K3kZoP3wQ8DNyTd60HvDFv+78R8U9JKwNXAL+OiDMqlf41RtK8ESNGjJg3b17rxGZmBsDIkSPp6+vri4iR9fs6rpHlAPUu4MvA7sAH8q77gZ8Cp0TE3Jz2BdIzNTMzs57wFFVLmWtkZmada1Yj8xRVZmY2pDmQmZnZkFZpGRdJWwCHA+8kTUWluiQREev2s2xmZmYtVZk0eAJwAfAScBdpcmAzM7MBUaVGdiwwC9ghIuZ0uTxmZmYdqfKM7G3ADxzEzMxsMKhSI3sIWKnbBTFr5Izpdw3Ieb+449gBOa+ZdaZKjey/gAMkLd/twpiZmXWqSo3sb8C+wF8lfR+4jzQx8GIi4g/9LJuZmVlLVQLZ1YWfzyOtAl2kvM01NjMz67kqgeyQrpfCzMysoiqTBk/pRUHMzMyq8BRVZmY2pFWaogpA0makJVxWZ8mAGBHxjf4UzMzMrB1VpqgaBlwCfJhFHTtqcy1GYZsDmZmZ9VyVpsXjSEHsJNKimQImArsA1wM3Aht2q4BmZmbNVAlkHwF+GRHHAbfnbf+KiN8DOwArAgd3p3hmZmbNVQlkbwGuyz/XBkKvCBARLwM/A/6t/0UzMzNrrUoge4ZFz9aeAV4F3lTY3wes1c9ymZmZtaVKIJsNjAWIiFeAO0jNjUgSMB54sFsFNDMza6ZKILsK2LcwafA5wM6SZgN3k56Tnd+l8pmZmTVVJZB9i0W9FYmIs4GjSE2KTwFfBU6tUhhJq0o6U9LDkuZLuknSnm0eu66kSyX1SXpG0jRJS/SelHSSpCskPSYpJB3fJM/3Srpa0nOSnpJ0kaQ3V7k2MzPrjY4DWUQ8GxGzcseO2rbvRMR7IuJ9EXFKRNRPJNyuqcABwCRgN2AmMFXSrs0OkjSG1PV/HdJQgP2BUcB1ktauS/7vwGrApS3y3ACYQQrYHwE+BWwKzJC0aicXZWZmvVN5Zo9uy8FqB2B8REzN264F3gGcDkxrcvhRpBlGNqutXC3pBtISM8cChxXSrhYRr0oaSQpOjZxA6syyR0Q8l/O8nfRM8PPAKR1fpJmZdd1gmmtxH1Lz5GW1DblmNwUYV9ZMWHfs9FoQy8fOBS4ndT6hsP3VVgWR9Dpgd+BXtSCWj70T+DNpPTYzMxsEBlMg2xiYWRJobi3sX0KeMmtdFg3Orj92TG567MQ7gGFN8iwtSy7PvGYvYESHZTEzsyYGTdMiMBq4q2T7k4X9ZVYnPcd6smRf8djHOixL8fj6PIdJGhYR8zvI04aYM6aX/TkuHV/cceyAndtsqBlMgQyWXG263X39PbZreUbEyGYZulZmZtZdLZsWJf1I0gcKv28taY0elGUu5bWuUfm9rHYEqct/VDy2WVlokuf8iHihwzzNzKwH2nlGdjDpGVTNtcCOPSjLHcAGkurLtEl+L3teRW7eu5fy51abAI9HRCfNiuT85jfJs7QsZma29LUTyJ4A1iz8rkYJ+2kqMBLYo277BGBWRMxsceyOkhbO8ShpVM7rkk4LEhEvAb8lzWCySiHPscAWVfI0M7PeaOcZ2Z+ASZLeSmrGAxgvab0mx1RZIXoaqbZ3vqTRpDFgE4Etgb1qiSTNALaJiGJAPQ04CJgm6QTgZdKg6peBk4snkbQNsAZQC1AbSvpIrQwR8Xz++evAX4FfSzoNGE5ag+1+4PsdXpuZmfVIO4HsP0hjuY5g0erP46kbn1Wn4xWiIyIk7U0KPCeTamczSQOkL29x7KOStiIFtAtJNc3rga0j4oG65CcA2xR+3y+/AN5OClRExExJ25EGPl8MvARcCRwZEc90cm1mZtY7LQNZRNwPbCNpRdLyLPeTgttlTQ6rJCKeBg7Pr0Zptm2w/W4KNbdOj2+Q9kZg+3bTm5nZ0td29/uIeBF4QNIU4C8R8c/eFcvMzKw9HY8ji4hDelEQMzOzKipNUSVpuKQTJN0q6dn8ulXS8ZKGd7uQZmZmjXRcI8vd2q8HNiB1zb8l7xoLHAfsJ2mriOh0ELKZmVnHqtTITgTGkTpkvDEitoqIrYA3kZY3WR84vmslNDMza6JKINsTOC8izo6IV2obI+KViPgB8CNg724V0MzMrJkqgWxNFjUnlrmZxWcCMTMz65kqgexRYNMm+zfNaczMzHquSiC7HDhU0meKE/xKWk7Sp4FPAL/uVgHNzMyaqbIe2XGk2e/PBk6QNCtvX580h+E9pHkKzczMeq7jGllEzAU2A75FWrfrffn1BPBN4H05jZmZWc9VWiE6z4l4bH6ZmZkNmEoze5iZmQ0WDmRmZjakOZCZmdmQ5kBmZmZDmgOZmZkNaZV6Ldqy6Yzpdw10EczMltBxjUzSXZKOkbRWLwpkZmbWiSpNiy+RBj4/IOlSSbsXp6oyMzNbmjpuWoyIjSRtDhwKfBTYA3hE0mTgRxExu7tFNFv2DFQz7hd3HDsg5zXrj0o1qYj4c0R8Cngj8EngPuArwF2SrpH0cUkrdbGcZmZmpfrVJBgRz0fEBRGxJWnV6IuAbYELgTmSzpD01v4X08zMrFy/n21JWl7SPsB3gI8BAVwL/Bn4AvAPSXv19zxmZmZlKgcySeMkfRv4F3AxaUb804CxEbFDROxGqqXNAk5tM89VJZ0p6WFJ8yXdJGnPNo9dN3c+6ZP0jKRpkjZskPaI3PtygaTZko4u67AiaV9Jf5L0VH7dIOmj7ZTHzMyWjird7z8h6Y/AHcCXgL8D+wFrR8SXi509IuIe4Exg3TaznwocAEwCdgNmAlMl7dqiTGOA64F1gInA/sAo4DpJa9elnQScQWoG3Qk4HzgJOLku3UTgV8Ac4OP59S/g55I+0eb1mJlZj1UZEH0e8AhpPbJzI+L+Fulnkp6ZNZWD1Q7A+IiYmrddC7wDOB2Y1uTwo4DVgc0iYk4+9gZSJ5RjgcPyttH597Mi4rh87AxJw4GjJZ0VEQ/l7YcA/wQ+GhGv5uN/D9wLTAB+1OqazMys96o0Le4LvCUijm0jiBERf42IQ9rIdx+gD7iscGwAU4BxjZoJC8dOrwWxfOxc4HJgfCHdzsDKOc+iyaSgXmzGfAl4thbEcp6vAs8CC9q4HjMzWwqqBLI9SM/DSkl6v6QqtZWNgZnFwJHdWthfdr5hpKbL20t23wqMyU2PtTyC1Cy6UETcDcyvO8dZwAaSjpX0BklrSDoWWJ/UNGlmZoNAlUB2MM2feb2d9JyqU6OBJ0u2P1nYX2Z1QG0eOxp4PiLKalRPFc8REZeRamhHAY8Dj5HGyu0XEb9rdBGS5jV7ASMaHWtmZp3rxdRSw0nNclVExX2dHNtWOkk7Aj8Ffg58GNiF1BnlZ5J2a1EWMzNbStrq7JEHNa9T2DRO0tYlSUeROlbcU6EscymvdY3K72U1Lkg1qWjz2LnAcEkrldTKVq+lkyTSc7RrIuKzhTS/y70g/xP4bVlhImJkg3KS83atzMysi9rttXgI8HVSwAhSz79jS9IJeDWn79QdwL6Slqt7TrZJfi97BkZEzJd0L+XP0DYBHo+IxwrnELARcPPCQkvrAcMK51iTNP3WTSV53gRsK2nliHihrSszM7OeaTeQXQrcTwoCPwJ+CNxQlyZIPfpujIgHK5RlKmki4j0o9FwkdXWfFREzWxx7uKS1IuIRAEmjcl4/K6S7gtTj8CAKgYz0TO9lUi9HSLW8F4D3l5xrc2Cug5iZ2eDQViCLiL+TBj4j6W3AxRFRWkPqh2mkqa3Oz+O97iMFmC2BhVNcSZoBbBMRKhx7Gik4TZN0AikoTcrvCwc6R8RcSd8EviapL59vC+AY4Lu1ABwRCyT9F/Afks4jDYxenhRUt8x5m5nZIFBlGZcTelGQiAhJe5MCz8nASNJg6vERcXmLYx+VtBUpoF1I6sRyPbB1RDxQl/xE0ni1z5N6Ic4hNZueUpfuKOBO4NPAR0hNpneRAuZPKl6mmZl1mdKY4yYJpAn5xwtzsJnQ9IAsIn7c38K9FkmaN2LEiBHz5s0b6KJ0bKDWyLKlx+uR2WA1cuRI+vr6+so61LVTI5tMev51EfBi4Xc1PoQAHMjMzKzn2glk2wFExIvF383MzAaDloEsIq5r9ruZmdlA6sXMHmZmZktNlWVcliBpBVIX+VHA5bWxXGZmZr1WZWHNUyXdWPhdwFXAL4BzgNsktbuQppmZWb9UaVrcmTRGq2YPYGvg26RVlAG+3M9ymZmZtaVK0+JbgLsLv+8B3BcRXwaQtBFwQBfKZmZm1lKVGtmKwCuF37cjNS3W3EuacNfMzKznqgSyB0kT59ZqX+8Ail3yx5AmDzYzM+u5Kk2LF5Em3R1DWg7ladKEvzWbArO7UDYzM7OWqtTIvkmapmoL0lRUEyJiHoCkEcCewNXdKqCZmVkzVWa/X0BaN+zQkt3PkJ6PPd/PcpmZmbWlKwOia/LKzn3dzNPMzKyZSoFM0hbA4cA7gdEsORN+RIQHRZuZWc91HMjyemQXAC+RFpqsX7jSzMxsqalSIzsWmAXsEBFzulweMzOzjlTptfg24AcOYmZmNhhUCWQPASt1uyBmZmZVVAlk/wUcIGn5bhfGzMysU1Wekf0N2Bf4q6TvA/ex+NyLAETEH/pZNjMzs5aqBLLirB3nkWb3KFLe5hqbmZn1XJVAdkjXS2FmZlZRlSmqpvSiIGZmZlVU6eyxkKSVJL1Z0ordKIykVSWdKelhSfMl3SRpzzaPXVfSpZL6JD0jaZqkDRukPULSXZIWSJot6WhJS9wLJZ+W9DdJz0uaJ+nPkj7Y32s1M7PuqBTIJL1H0jWkSYIfALbM28dIulrSDhXLM5W0uvQkYDdgJjBV0q4tyjMGuB5YB5gI7A+MAq6TtHZd2knAGaTlaHYCzgdOAk4uyfo84FTgYmDXXLZpwPBKV2dmZl1XZYqqd5OCxhPAjyk8M4uIxyQNIwWTq8pzaJjvrsAOwPiImJq3XUtauPN0Fl/zrN5RwOrAZrWB2pJuIPWoPBY4LG8bnX8/KyKOy8fOkDQcOFrSWRHxUE67L3AwsGVE3FA41287uS4zM+utKjWyE4E5pEU1v8ySEwZfDby/Qr77kGbOv6y2ISICmAKMa9RMWDh2enG2kYiYC1wOjC+k2xlYOedZNJkU1IvNmF8A/lAXxMzMbJCpEsi2As6NiGdZsus9pKbGN1XId2NgZl4KpujWwv4l5BrgusDtJbtvBcbkpsdaHgHcUUwUEXcD82vnkPQ6YHPgNkknS3pU0suS7pA0sfNLMzOzXqnS/X5lmq85tlrFsowmzaZf78nC/jKrk2qFT5bsKx77WH5/Pi8OWu+pwjlGk6bhmkiakutwYB5pMdHJklaMiHPLCiNpXoNy1oxosd/MzDpQJZDNBt7bZP/2pE4aVZTV8NrZ18mx7aSr1VRXBnaNiH8CSLqK9MzuOKA0kJmZ2dJVpWnxp8BBdT0TA0DSkaTnUBdWyHcu5bWuUfm9rMYFqSYVbR47FxguqWzS49UL6Wp53lkLYrDwmd3vgLULzZWLiYiRzV54BW0zs66qUiM7DdgR+D1wJ+kD/wxJawBrAdOBsyvkewewr6Tl6p6TbZLfy56BERHzJd1L+TO0TYDHI+KxwjlE6qhycy2RpPWAYbVz5DzvaVDOWueW+md5ZmY2ADqukUXEi6RAdhSpg8QLwFhSd/yjgd1LOmy0YyowEtijbvsEYFZENGuunArsKGmt2gZJo3JelxTSXQEsAA6qO34i8DKpl2PNJcAGktYp5ClgF+DeiHii9SWZmVmvVamREREvkwYVn9HFskwDrgXOz+O97iMFmC2BvWqJJM0AtomIYrf/00jBaZqkE0hBaVJ+XzjQOSLmSvom8DVJffl8WwDHAN+NiAcLeX6bNAD6dznPWmeP9wL/1sXrNjOzfqgUyHohIkLS3qTAczKpdjaTNED68hbHPippK1JAu5BU07we2DoiHqhLfiLpOdXnga+QxsR9HTilLs+5Oc9vk5pKhwG3AftExKX9uVYzM+sepf4LtrRImjdixIgR8+a16qU/+JwxvWx0hL2WfHHHsQNdBLNSI0eOpK+vry93mltMyxpZnlOxUxERH6pwnJmZWUfaaVp8B0uOvRoOvCH/PI/Uk6820PcJ4NmulM7MzKyFlr0WI2KdiHh77QV8iNRb8XvAmyJiVESsTpqW6kzg+ZzGzMys56oMiD4D+FNEfDEiHqltjIhHIuI/gD/T3d6MZmZmDVUJZNsC1zXZPwPYrkphzMzMOlUlkAWwQZP9G9F6XkQzM7OuqBLIrgQOkzQhz3QBpFkv8hInn8lpzMzMeq7KgOgvAe8DLgC+JeluUg1sLLAm8GBOY2Zm1nNV5lp8CHg3aSaMp0irQX8g/3wK8O6cxszMrOeqzrXYB3w1v8zMzAZMlWdkZmZmg4YDmZmZDWkOZGZmNqQNmmVczGzgDeQKB55536pyjczMzIY0BzIzMxvSOg5kku6SdIyktXpRIDMzs05UqZG9BHwTeEDSpZJ2l+SanZmZDYgqM3tsBHwQmEKa5f4y4EFJJ0lat8vlMzMza6pSTSoi/hwRnwLeCHwSuA/4CnCXpGskfVzSSl0sp5mZWal+NQlGxPMRcUFEbAmMAy4irVd2ITBH0hmS3tr/YpqZmZXr97MtSctL2gf4DvAx0kz415JWiv4C8A9Je/X3PGZmZmUqBzJJ4yR9G/gXcDGwGXAaMDYidoiI3Ui1tFnAqd0orJmZWb2OZ/aQ9AngUGDzvOkq4IfAZRHxcjFtRNwj6UzgvP4W1MzMrEyVKarOAx4BvgWcGxH3t0g/k/TMzMzMrOuqNC3uC7wlIo5tI4gREX+NiEPayVjSqpLOlPSwpPmSbpK0Z5vHrpvHtfVJekbSNEkbNkh7RB7YvUDSbElHNxsLp+QaSSHpu+2Ux8zMlo4qgWwP0vOwUpLeL+lHFcszFTgAmATsRqrNTZW0a7ODJI0BrgfWASYC+wOjgOskrV2XdhJwBqmH5U7A+cBJwMlNTvEp0vM+MzMbZKoEsoOBZgOf304KJh3JwWoH4JMRcX5EXJPzuQE4vcXhRwGrA7tGxKUR8RtSIFwJOLZwjtH597Mi4riImBERJ5M6oxxZH/TyMW/O+7/Q6TWZmVnv9WIZl+Gkaaw6tQ/QR5opBICICElTgB9K2jAiZjY5dnpEzCkcO1fS5cB44LC8eWdgZdKsJEWTga8CewJn1+37AfCHiLhYUoXL6q6BXGbDzGwwaiuQ5UHN6xQ2jZO0dUnSUaSgcU+FsmwMzIyIV+u231rcX1K2YaQa4i9L8rwV+LikMRHxWM4jgDuKiSLibknz8/5i3vuTpuEqfdZWRtK8FklGtJuXmZm11m6N7BDg66QgEKTmuWNL0gl4Nafv1GigrLrxZGF/mdXzeZ8s2Vc89rH8/nxELChJ+1TxHJLeAHwPODYiHmxZejMzGxDtBrJLgftJAeNHpHFjN9SlCeBZ4MZ+fPBHxX2dHNtuujNJc0ie1eK8i2cQMbLZ/lxjc63MzKxL2gpkEfF34O8Akt4GXBwRt3e5LHMpr3WNyu9lNS5INalo89i5wHBJK5XUylavpZO0I2m6re2B1eqeja0kaSTwbP0AcDMzW/qqLONyQg+CGKTnVhuUjOfaJL+XnjMi5gP3Uvd8q3Ds4/n5WO0cAjYqJpK0HjCscI6NSPdmBilQ1l4An80/79DORZmZWW+1rJFJmpB/vDD3IpzQ9IAsIn7cYVmmkqa+2oNCz0VgAjCrSY/F2rGHS1orIh7J5R6V8/pZId0VwALgIODmwvaJwMvA5fn3XwH/W3Kea0nzSp7Fok4oZmY2gNppWpxMarq7CHix8HuzvugBdBrIppECxfl5vNd9pACzJbBw9nxJM4BtIqJ4/tNIwWmapBNIQWlSfl840Dl3yf8m8DVJffl8WwDHAN+tPduLiIeAh+oLmJsYH4qIGR1em5mZ9Ug7gWw7gIh4sfh7t+Xa3t6kwHMyMJLU3X58RFze4thHJW1FCmgXkpoFrwe2jogH6pKfSBqv9nnSYqBzSD0yT+ni5ZiZ2VKiiFadAa2bJM0bMWLEiHnzWg03K+cB0fZa9cUdxw50EWwQGzlyJH19fX1lPcP7vbCmmZnZQOqks0dHKnT2MDMz61gnnT06mWiwSmcPMzOzjrXd2cPMzGwwahnIIuK6pVEQMzOzKtzZw8zMhrTBNLOHmZlZxwbTzB5mZmYdGzQze5iZmVXRcWcPd/4wM7PBpN2FNZcgaSVgW+AdedO9wHUR8UIXymVmZtaWSoEsd/j4DmkxytqzsgDmSToyIiZ3p3hmtqwYqHlEPcfj0NdxIJP0MVKHjwdIs83PJAWzDUmLTp4vaX5E/LyL5TQzMytVpUb2VeBOYPOIeLqw/TJJZwN/AY4FHMjMzKznqgyIXh+4oC6IARARfcAFwDv7WzAzM7N2VAlkj9B8DNmrwKPVimNmZtaZKoFsMnCwpFXrd0haDfgEqVZmZmbWc+1MUbV13aY/ALsDt+VnYneSeixuCBwGPAFc3+VympmZlWqns8cMUqAqqjUtnlLYV9v2NmA6sHx/C2dmZtZKO4HskJ6XwsxsgAzU+DXwGLZuaWeKqilLoyBmZmZVeD0yMzMb0voz1+JmwAdI01TVB8SIiG/0p2BmZmbtqDJF1TDgEuDDpA4exbXJorDNgczMzHquStPicaQgdhJpbTIBE4FdSN3ubyR1xe+YpFUlnSnpYUnzJd0kac82j11X0qWS+iQ9I2mapNJySDpC0l2SFkiaLeloScvVpfmkpF9L+mcuy925bGtUuTYzM+uNKoHsI8AvI+I44Pa87V8R8XtgB2BF4OCK5ZkKHABMAnYjTUg8VdKuzQ6SNIYURNchBdX9gVHAdZLWrks7CTiDtOL1TsD5pKB8cl22JwBPA18BdibN9v9R4EZJIyten5mZdVmVZ2RvIX2oA7yS31cEiIiXJf2MNDD6K51kmoPVDsD4iJiat11LWu/sdGBak8OPIj2r2ywi5uRjbwDuI01gfFjeNjr/flYOxAAzJA0HjpZ0VkQ8lLdvGhGPFc5xnaSZpHF1BwH/2cn1mZlZb1SpkT3DogD4DGluxTcV9vcBa1XId5987GW1DRERwBRgXKNmwsKx02tBLB87F7gcGF9ItzOwcs6zaDLpmhY2Y9YFsZob8/vaJfvMzGwAVAlks4GxABHxCnAHqbkRSSIFjgcr5LsxMDMiXq3bfmth/xJy55N1WdTMWX/smNz0WBEQukUAABolSURBVMsjcpkXioi7gfmNzlGwfX4vO5eZmQ2AKoHsKmBfSbUpqM4BdpY0G7ib1Dx4foV8RwNPlmx/srC/TG2V6naOHQ08HxELStI+1eQcSBoFnEm6xl80STev2QsY0ehYMzPrXJVnZN8CLiR3uY+IsyWtDBxIemZ2LnBqxfLUz+nY7r5Oju34HJJWAS4ldSDZukEgNDOzAdBxIIuIZ4FZddu+w6IOIFXNpbxGNCq/l9W4INWkos1j5wLDJa1UEoxWLztHbrr8NbApsFNE3FqfpigimvZodK3MzKy7BtMUVXcAG9SP5wI2ye+lz6UiYj5wL+XPtzYBHi903LiDVJPcqJhI0nrAsPpz5JrmZcAWwO4R8ae2r8bMzJaKyoFM0kqSdpJ0WH7tlD/4q5oKjAT2qNs+AZgVETNbHLujpIW9JfMzrT1Is5DUXAEsIHWfL5oIvEzq5Vg7fiVSc+JWwF4RcV1HV2NmZktFpbkWJU0gNSXWOlpAat6bJ+nIiJhcIdtpwLXA+Xm8132kALMlsFfh3DOAbSJChWNPIwWnaZJOIAWlSfl94UDniJgr6ZvA1yT15fNtARwDfDciir0tf0UaMH0i8KykzQv7Ho+I2RWu0czMuqzKXIsfI427eoAUQGaSgtmGwGdJgWh+RPy8k3wjIiTtTQo8J5NqZzNJA6Qvb3Hso5K2yuW5kFTTvJ7UMeOBuuQnksarfZ40aHsO8HXSIqFFu+f34/KraArVZy8xM7MuUhpz3MEB0t+B1wGbR8TTdftGAH8BXoyI/9O1Ur6GSJo3YsSIEfPmzat0/EAuAmhm3eWFNds3cuRI+vr6+so61FV5RrY+cEF9EAOIiD7gAuCdFfI1MzPrWJVA9giLnouVeRV4tFpxzMzMOlMlkE0GDpa0av0OSasBnyDVyszMzHquZWcPSVvXbfoDqSPEbZLOBu4k9VjckDTL/BOkjhZmZmY9106vxRksOXVTrWnxlMK+2ra3AdOB5TEzM+uxdgLZIT0vhZmZWUUtA1lE1K/dZWZmNmgMprkWzczMOlYpkEkaLukESbdKeja/bpV0vKTh3S6kmZlZI1WmqBpF6pW4AamH4i1511jSVE77SdoqIhotu2JmZgzcTD2vtRlFqtTITgTGAYcDb4yIrSJiK+BNpPkL1weO71oJzczMmqgSyPYEzouIsyPildrGiHglIn4A/AjYu1sFNDMza6ZKIFuTRc2JZW7OaczMzHquSiB7FNi0yf5N8VyLZma2lFQJZJcDh0r6jKSFx0taTtKnSXMt/rpbBTQzM2umygrRxwE7AmcDJ0ialbevD6wB3ENaqNLMzKznOq6RRcRcYDPgW8Bc4H359QTwTeB9OY2ZmVnPdVQjkzQM2A+YFRHHAsf2pFRmZmZt6rRGtgA4j+adPczMzJaajgJZRLwKPACs1pvimJmZdaZKZ48pwEGSvhcRC7pdIDMz662BmhoLejM9VpVA9idgPPC/eYXou4Hn6xNFxB/6WTYzM7OWqgSy6YWfv0f56tGBV4g2M7OloEog+wRLBi8zM7MB0XEgi4jJPSgHAJJWBU4mdfEfCdwBnBgRLWcKkbQucDqwHakTy/XAURExsyTtEaTZ+98GPAScA5yWO7NUytPMzAZGR70WJa0h6QP5A74XpgIHAJOA3YCZwFRJu7Yo1xhSkFkHmAjsD4wCrpO0dl3aScAZwEXATsD5wEmkAFopTzMzGzht1cjynIpnA58kPQND0g3APhHxeDcKkoPVDsD4iJiat10LvINUK5rW5PCjgNWBzSJiTqF895EGbR+Wt43Ov58VEcflY2fkVa2PlnRWRDzUSZ5mZjaw2q2RHQ58GngEuAS4DfggqUmuW/YB+oDLahsiIkjd/cdJ2rDFsdNrAScfO5c0wfH4QrqdgZVznkWTSUF9zwp5mpnZAGo3kE0A/gFsEBH7RcS7SU1ye0ga2aWybAzMrH9OBdxa2L+EPG3WusDtJbtvBcbkZsJaHkF69rZQRNwNzK+do8M8zcxsALXb2WN9UqeLZwrb/hM4FBgL/LULZRkNlI3Se7Kwv8zqpObOJ0v2FY99LL8/32Ag91OFc3SS52IkzWtQzpoRfX19jBxZLf4veLk+zpuZDR0nrFBl9TDo6+uDBrNKtRvIhgNz6rbNKezrlmbd+lt1+W/32E7O0Z/yNBN9fX1P9+P4oWhEfu8b0FK8Nvhedo/vZfe0dS9fqJ7/akDpN/lOut83+pBXlRKVmEt5rWtUfi+rHUGqSUWbx84FhktaqaRWtnohXSd5LiYiutXU+ppSq6n6/vSf72X3+F52z0Dey04C2a6S1ir8vgrpw34/Se+uSxsRcUaHZbkD2FfScnXPyTbJ72XPq4iI+ZLupfwZ2ibA4xFRawK8gxR4NwJuriWStB4wrHaODvM0M7MB1Ekg+3h+1ftMybYgjdXqxFTSM7c9KPRcJHU0mdViEPJU4HBJa0XEIwCSRuW8flZIdwVpKZqDKAQy0jixl0k9EjvN08zMBlC7gWy7npYimQZcC5yfx3vdRwowWwJ71RJJmgFsExHFJs3TSMFpmqQTSEFpUn5fONA5IuZK+ibwNUl9+XxbAMcA342IBzvN08zMBlZbgSwirut1QSIiJO1NChInk6aomkkaIH15i2MflbQVKfhcyKLppLaOiAfqkp9Iehj5eeArpE4rXwdO6UeeZmY2QJTGHJv1lh+qd4/vZff4XnbPQN7Llh36JX2oauaSdqh6rJmZWTvaGZn2O0nXSNpdUss1xiS9TtI+kq6j+fyIZmZm/dbOM7JNge8AvwaekDSdNJPHbNJYKpHGVr0T2BzYnjQm60qgvlu+LaPcdNM9vpfd43vZPQN5L9t+RiZpC+BzpB6Eq1K+MvTTpEmFfxARN3axnGZmZqU67uyRmxffC2wIrEEKaI+TBhPfUjLpr5mZWc+416KZmQ1p1aYhtmWapLUlfU/S/0h6VlJI2rZB2o9L+rukFyQ9JOlbklYuSbempCmSnpD0nKTrJX2w5xczwNq9l5Luz/vqX98qSbus3ssPSZosaZak5/Pf2yWSNilJu6OkP0uaL+kxSeeULUklaVVJZ0p6OKe9SdKe9elea9q9l5JmNPi7vKgkz57dSwcyq2I9YH/gWeDqRokkHQj8BPgjsAtpoPvnSQuZFtOtnPPZBvgCaVHTZ4CrJW3a/eIPKm3dy+wPpJloiq/vFxMs4/fys8BbSdPj7QJ8Kf9+o6TNa4nyF4VpwIOkKeeOIi2q+1tJ9Z+JU4EDSLP67EaapGGq0or2r2Vt3cvsbpb8u5xUkmfv7mVE+OVXRy9gucLPe5Oek25bl2Z54GHgsrrtn8rpP1DY9rm87T2FbSsB9wJXDPT1DvS9zPvuBy5tI79l+V6OKdk2krSaxcWFbX8Fbqm79zvm+/axwrZd87Z9CtsE/A/wj4G+3kFyL2cA/9tGfj29l66RWceivQ49mwNrAVPqtv8EeAnYt7BtH+C2iFg4kXOkZXZ+Buwo6fX9K/Hg1ea97MSyfC+XWJEiIuaRagxrA0h6M/A+4MLivY+I6cC/WPLvso/CJOaRPoGnAOMkbdiDyxgU2rmXHerpvXQgs16pLYGz2PI7EfE8aQzixnVpy5bpuZVUs9ugFwUcgrbPz9FelHSbpMMk1a8H6HtZIGkNFr8npX+X2W0s+Xc5s+TLxq11eS0TSu5lzfqSnpL0sqS7JU2S9Lq6ND29l50s42LWidqipGULkD7J4ouWjm6SDsoXOF3W/Aa4idREOBo4EDgbGAt8sZDO9zLLQf6HpC/sp+XNrf4u31P4fTRwV4N0xbxe8xrcS0gTqV8E3EkaX7w3aWL295JqYTU9vZcOZNZrjcZ3NFpxvJM8lhkRcXjdpqmSfgIcIem7EfHPYvJmWXW/dIPWt0kfrIdExD/q9vnvsjOl9zIivlaX7jeSHgW+KmnLiPifwr6e3ctKTYuSxub5FD8j6dP553f2pyD2mjM3v5d90xrF4t+I5zZJB+Xfni09X1gOeH9hm+8lIOkk4Ejg3yNicmGX/y471OReNlJ7Lr5FYVtP72XbNTJJG5C6ZO4HrFnbnN8jp3kU+AVwTsk3IFu23JHfN6bQpCBpFWBdFl+N+w7K28g3AV4hNVvYkmpfRIvPHZb5eynpROCrwNERcWbd7uLf5ZV1+zYB/lSXdl9Jy9U926mNpSp7zvaa0uJeNtLo77Jn97KdZVzWlfSrfKJDgb8DJwATSF0qd8s/n5j3fRK4XdIvJb2jP4WzIe3PwCOkVbaL9gdeR5qTs2YqsImkhZNMS1oxp70qIp7ucVmHqgmkD4vivKbL9L2U9HXga8DXIuLb9fsj4iHSs8YDimPGlJarejNL/l2OJI01K5oAzIqImV0u/qDS6l42MSG//7mwrbf3so3+/wtI//AHAcPbSD8cmAj8DXhhoMdD+NWbF/CR/DqFVCP/ev59l0KaiXnfWcC2wGGkiaV/WZfXyqTBkfcCHyON6fkNMB9470Bf60DfS1IQuij/H9yO1EV8ak57qu/lwms/Mt+Ty0nDP4qvTQvptgdeBn4OfCjf1zmkD97lC+kEXAM8AXwi3/vJpC8Pewz09Q70vQS2An6b782HSIPKz8/35xd1+fX0XrZzQXv142ZUPtavwf3Kf+Rlr/vr0h1I6ta8gDRO51RgWEl+awEXktrKnycNlNxyoK9zMNzL/OFxFWmA+YukmTr+BExskN8yeS9Jg3Pb/bvcGfgL8AJp0vNzgdVL8lyN9EXskZz2ZmDvgb7WwXAvSbPS/BZ4KN+b50kDzf+DwheCpXEvPWmwmZkNaR4QbWZmQ1rHgUzSWZJultQn6SVJcyT9JnfFX60XhTQzM2ukysKar5La6efk99EsGh/wFPD5iFhiCn8zM7NeqNK0+C5glYh4R0SMi4g1SA/9vkR62Pffy8ASB2ZmNkh0tbNHnln7emBBRHygaxmbmZk10NXOHhHxDHAei0Zrm5mZ9VS/Jg2WtALwb8As0lxaY0hTWD3a/6KZmZm11t/Z71cCfsyimYtfIA2O+0I/8zUzM2tLv5oWI+I5YH3gM8DvSYFsSkT8pgtlM1tqJG2cFwbccaDL0g2StpUUkg4e6LL019K6ll6dR9LeeTFUrxDSI1XGkS02Viwi7o6I8yJiV2AH4DBJx3WrgGZLyXeAP0Za8n4JklaX9EL+oDtwKZfNBoikd0s6XtI6VfOIiEtJ07Sd0q1y2eKq1MhmSTq4OHN0TUTcQvrH+nS/S2a2lEjagjS57neaJDsAWBG4j7QKhL32/AEYRpqnsubdpEmc1+ln3t8D9pG0UT/zsRJVAtn9wI+AuyQdJWn92o7c+WNbYPWulM5s6fgcqbPStCZpDgWuBb4LbCNp3aVRsMFE0vJ5PbnXpIh4NSJeiIhXepD9JaRxtp/tQd7LvCqB7IOk5ThWIc1kPlPS85LuJc3ssQ9pLJlZT0haTtINkh6VNLJBmv+TmwHPbpHXCqQl3KdHxEsN0ryH9M18CvAT4CXgkAZpD87n3T5/0ZstaYGkuyRNLEm/jqSLJT2dp327TNLbJd0vaUYh3fE533VK8lgsbZNrfb2k/yfpL5KeyOW6R9K36gNU4Tp2kPQ1SbNJz8A/2iDvXXL6Ixrsv0HS45JeV9i2kqSvSrojN9vOk3S5pE1bXUs+/g2Svi/pwfwM6sH8+xIrEUtaUdLRkv43f171SbpJ0uGFNIs9I5N0PHBB3n1t3heSJksan3/+ZIOy3ZHvrQAi4lnS5+J+7VybdabjXouRRlCfI2kyqev9PsCmwNqkP/TLgc93sYxm9T5FWtrkkxExr0GaO0hTqH2wRV7vBVYF/tokzaHAc8DFEfGcpN8CEyUdF4uvdlt0MqmZ6hzSEjaHAZMl3RMRfwTIH7jXk1Zc/y/gH6Q1nq4lrevXbW8mLXx7MfBT0ppc2wBHk/4P71RyzGmkhVDPJa0lN6tB3leSlpmZACy2knDu5LA5cGbty0IOaL8j/ftcSFreYwTp3/aPkraOiJsaXYikEaSlbNYjtRDdnK/hMGB7Se/P41prC4v+ntRadCXw36TPqk2A8fncZS4B3kh6VHIy6d8HYDZpMdNHSH8b59WVbXNgQ+DYWHzGiRuAnSSNi4jX/ErdS9VAr3vjl1+dvEgLRz5MWjhyhRZp7wOeaZHmENLwkT2bnO9JYHJh2175mF1K0h+c990CrFjY/mZSQPtZYdupOe0BdXnUts8obDs+b1un5Jz3F9Pmbdvm9AcXtq0IvK7k+G/ktO8vuY5ZpCnp2vm3+XY+ZsMG+b+nsO2LedtOdWlXAx6ou/ayazkpb/tc3fGfz9u/Udh2dN52ckmZl2txntp92Lbk2JMbXO+5pC8Jb6rbfmBOv+9A/z96rb28jIsNNR8jLRx5TkS8XNsoadVis1X2IqkJvJk18vuTDfaPJz3znVLY9lvgMdJKt42cHREv1n6JiH8BdwHFLth7kILyz+qOPa1FmSuJiBdjUY1ohdwT8w2kRTsByqaV+0FEPN/mKWr3qLbUPblp7UDg9oi4uZD2QOBO4G+5ifANuSwrAtOBLSUNa3KufUgLYv6wbvs5pFWI9ylsO4D02OPE+kyicY26HeeSAtPCzj+ShpP+Rq+IiDl16efm9zH9OKeVaBnIJH2oauaSdqh6rFkDtQ+oX9c2SBpLWjV5/7q0o0kfds3Umn7UYP+hOY+HJK0naT1SD7bpwJ75w7fMvSXb5rJopQiAtwP31H+YRsRjQKMm036R9DlJt5Jqh0+Srm1G3l3WSeuudvOOiNtJNdEDtKhX89ak+zWlLvkGwLh8/vrXJ4DlgUb3FtK9m1X8MpPL8DKpFvmOwuZ3AndGxAvtXks7IuI+0peAgwpfoj4KvJ665sas9jfm1Yy7rJ1nZL+TdD2pa/IV0aJHT/4H3Z203PUWpG9YZt3yHtKHcPEZw9b5feHzm9yrcDRwRYv8aoFuVP0OSW8HtiN9ADX6QD+Q1JOxXqP/J40CZivNPvzaetYt6UvA6aTnRGeyaCmmNwOTKf9i225trGYK6X5sT/qQn0C6Fz+pLw5pbNWXmuTV6ktIJ3oVPH4I/BLYk/Ts8VDSs7PflqSt/Y1187qM9v4DbEoKYr8GnpA0nfRgfDbpG51I/0C1B7rbk77ZXUnq6WXWTW8E/hX5oUM2Pr/PLWzbI7+3mmXm9vxeNuvCIaS/709RXkP6f6QPrrJA1o77gfUkLVeslUkaA9T3xqw1fY7Kx9XSrky6J/e0cb6D8rG71J1v5wplb+SnpGdlEyT9EfgIqUfow3Xp7iY1615TsXnvXmB9SSvUNTGvAIxl8RrxXcAGklaKiAUdnqdVALyM1Mx8qKTbgf8LnFJfU8zWy++3l+yzfmjZtBgRt0fEh0n/QL8nfUCcQfoHvJ40iPBSUrv+h0kBb/OI2CUiZvaq4LbMegFYI/dEI3fV3jxvH5O3vR44khR8ftoiv1tIvfE2L27MTWMHA7dFmrnmV/Uv0rOtjSW9r+K1XE4KQvVNokeVpK3VCOub679I+8NoXiF9MC+sFeYP/i+3eXxLEfE4qRY8nvRsajWWbFaENEfrWjSokUlas8WpLiUFwvru75/K26cWtv2E9OV6Usl5WtWQn83vS9TYAfIzx8mkHp9fz5vPb5DX5sCjEdGo56dV1Hb3+4i4AbhB0vKkLssbkv5gglRVvh24pZ8PT81auYbUjPOzPHbqGFKLwUTgW5J+SaolvRk4KBp3zwcgIl6RdAmwV9039g8Db6HxhxKkpqTj8/lurHAtpwAfBy6Q9H5Sc+mWpC+NT7B4beCqvP/E3G3/vpx285y2Hb8Cvglcka95tXz+0vFz/TCF9G90OtBH+tJb73uk2VS+LWl70r/r08BbgQ+Rvphs1+Qcp5LGZH1faZzfLaTWo0NJTcyn1p1rD2BS/tJxZc5/I9Jcsc2e5d8IvAocK2l10jCM+yLiL4U05wL/H+kLyXURcXd9JpJWJQ2t+FGTc1lVA91t0i+/OnmRxiteCcwnPeM5gVTD2AP4J+kD6kZgjw7yfD913aJJzz0C2KTFsbNINb9h+feDadxdewZwf922t5PGKz1D+iC/LG97AphWl3YsaezV8/mcvyAF7Ptpr/v98sBXSM2QC/L9OpXU8SKA4wtpG15HG/dzRVIzbwDnNkm3AnBE/vd6Lr/uJtWgPtzsWvL2NYCzSStuvJTfvw+8oeRcKwPHksYXvpDv340Uuu83Oc9EYCbpeWJQGIpRSHN13ndQg2udmPdvPND/h16Lr0orROdq/0dJvZGeJQ1GnB7td9M1G1Qk/Q4YHhFbDYKyjCYFsnMiwlMaDQGSppE6t70pIuaX7P8b8M+IGL/EwdZvHc/sIWkr0px0q7B4D6y5kr4REWeWH2k2qB0J/F3ShyPiyqV1UknDSj74jsnvpTPx2+CSh2TsBHy/QRDbmzSLyL8t7bItKzqukUn6C+n52OGk6vRypHb6L5GaaC6KiI93uZxmr0n5Od8/gZtITX8fIg1f+ROwdfRmAlvrAkkfIDXLHpHfN4iI+we0UMuoKoHsOeD0iFhizTFJh5LGVfx7RDSav8zMMklHksZarUOam/Eh0jOzEyLPFWiDU55vdgKpq//REXHJwJZo2VUlkD0CTIqIspHrSPoJ8K6I2LgL5TMzM2uqylyL1wK7tti/zK3VZGZmA6NKIPsh8EFJ/95g/zqkbtFmZmY9V6Vp8VXSEgXLk2b2OA/4G4vWNjoDOCYimi5oaGZm1g0dd78nrS30LtI8ilvnVzEa3gLMk7QxacbpsjnHzMzMuqLSgOiFB6dl5jclBbV355/HkQJkkEbbzwJujYiD+l1aMzOzOv0KZKUZpslcN2ZRYNuUNC1L/WzeZmZm/db1QGZmZrY0Vem1aGZmNmg4kJmZ2ZDmQGZmZkOaA5mZmQ1pDmRmZjakOZCZmdmQ9v8DX1mzHbbeOUMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAFnCAYAAAAL7CZWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZxcVZn/8c9XNAlG6RBkVZE17DMgoMiPsMeVHREXhARcwAVFZ0TEAcRBNgfE8QeyBxgFt7C5gGELQdkEBtl3CIogCaRD2Eme+ePcIkWluvrWrVtd3VXf9+tVr+q699xbT3cl/fS595znKCIwMzPrpDd1OgAzMzMnIzMz6zgnIzMz6zgnIzMz6zgnIzMz6zgnIzMz6zgnIzMz6zgnIzMz6zgnIzMz67g3FzlI0gRgPWA5IICngTsj4oESYzMzsx6hvOWAJK0D7A/sASxf2Zw9V07yFPBL4NSIuKfEOM3MrIsNmowkrQ4cC+wKvAjMBK4HHgLmkBLSeGANYDNgIrAkMA04OCIeblfwZmbWHfIko5eBO4CTgGkR8fwg7ccCHwcOBNaLiDElxWpmZl0qTzLaOSIuLnTyFo41M7PekfuekZmZWbu0NLRb0mhJ75Q0qqyAzMys9xRKRpLeK+kq4DlgFrBFtn05SVdK2r7EGM3MrMs1nYwkbUgaUbc6cG71voj4J2kk3T6lRGdmZj2hSM/oSOAJ0qTXb7NorlHFlcD7WozLzMx6SJFkNBE4PSLms2iya7VZwEotRWVmZj2lSDIaA/Q32L9UwVjMzKxHFUlGDwEbN9i/LXB3sXDMzKwXFUlGPwc+WzNiLgAkfRP4MHBeCbGZmVmPaHrSazan6HJgS+BeYG1SuaBlgRWA6cBHI2JhuaGamVm3arpnFBGvAJOAfyMVTn0JmADMBr4F7OBEZGZmzXA5IDMz6ziv9GpmZh036EqvWdmfZkVEbFfgODMz60F5lh1fjcUnt44F3pF9PZdUhaEvez0bmF9KdGZm1hMGvUwXEatExKqVB7AdaeDCScBKETE+IpYmVV34MfBC1sbMzCyXIkO7LwZeiIhPDbD/AmBMROxSQnxmZtYDigxg2BqY0WD/NcA2RYIxM7PeVCQZBbBOg/3rUb+AqpmZWV1FktEfgQMk7S3p9eUjlOwDfDFrY2ZmlkuRe0bvIi2utzLwFPAAqSc0AVgeeBzYIiL+Vm6oZmbWrQpVYJDUBxwM7Ewa+g3wMHAxcFxEzC0tQjMz63ouB2RmZh3nckBmZtZxeSowvEHO8kAuB2RmZrk1nYyoXx7ozcCKpJ7WbOD5FuMyM7MeUto9I0mjgW8AU4CtIuIfpZzYzMy6XukDGCSdB7x5oHJBZmZmtdoxgOE64ENtOK+ZmXWpdiSjVYFRbTivmZl1qSKj6VYeYNd4YHvgQFKxVDMzs1yKlANayMCFUAXcC+wUEQ+2GJuZmfWIIkO7j2TxZBTAM8D9wBURsbDVwMzMrHe4HJCZmXWcywGZmVnHNZ2MJJ0l6f0N9r9P0lmthWVmZr2kSM9oMrB6g/2rAvsUisbMzHpSOy7TjQVebcN5zcysS+UaTZfNLVqlatPakras03Q8cADgYd1mZpZbrtF0kg4HDmfg+UWvNwUWAlMi4rzWwzMzs16QNxn9K7AhKdmcBZwGXF/TLID5wM0R8XjJcZqZWRfLdZkuIm4HbgeQ9B7gNxFxZzsDMzOz3uFJr2Zm1nGD9owk7Z19eV5ERNXrhiLi3JYiMzOznjFoz6iqMOqSEfFK1Ws1OCwiYonywjQzs26W557RNgAR8Ur1azMzs7L4npGZmXWcC6WamVnHNTOAoSkewGBmZnk1M4Ch0YCFWh7AYGZmueUewGBmZtYuHsBgZmYdl6scUD2SRgNbA6tlmx4GZkTESyXEZWZmPaRQzygb1HACsDSL7iUFMBf4ZkRMLStAMzPrfk0nI0l7AucDs4CfAneTEtK6wP7Au4BPR8Qvyg11eJD0GmlI/LxOx2JmNoIsBSyMiLpX5Ioko9uBtwCbRcS8mn19wI3AKxHxL8XiHd6y0YXq6+vrdChmZiNGf38/pJHWdee3FrlntBbwH7WJKHuXfklnA0cUOO9IMa+vr69v7ty5nY7DzGzEGDduHP39/QNeUSpSgeFJGs85Wgg8VeC8ZmbWo4oko6nAZElvq90haSlgX+DsFuMyM7MeUuQy3UxgB+AOSScD95JG0q0LHADMBmZK2rL6oIi4tsVYzcysSxUZwLCwZlPlBKqzrbK9a8oDSZrre0ZmZs3J7hn1R8S4evuL9IymtBiTmZnZGzSdjCLinHYEYmZmvcvrGZmZWcc5GZmZWccVSkaSPi3pT5L+KWlBncdrZQdqZmbdq+l7RpK+C3yPNLH1z8CzZQfVzU6cfn/H3vugSRM69t5mZo0UGU33JeAa4MMR8Wq54ZiZWS8qcpluKeCXTkRmZlaWIsnoNuDdZQdiZma9q0gy+i6wv6T3lh2MmZn1piKTXmdI2g+4QdL1wKPAgsWbxX4lxGdmZj2gyGi695Mqd78ZmJg9agXgZGRmZrkUuUx3EvAqsDMwPiLeVOfRFUVRzcxsaBQZ2v0vwBERcWnZwZiZWW8q0jP6J/BK2YGYmVnvKpKMzgL2klSkV2VmZraYIgnlOtJKrzdkK70+wuKj6byyq5mZ5VYkGV1R9fUZvHFVV8hWdgU8iMHMzHLxSq9mZtZxXunVzMw6zovrmZlZxxUeESdpE+D9wNIsntQiIr7fSmBmZtY7ipQDWhKYBnyQRYMVlO2Oqm1ORmZmlkuRy3SHkRLRUcA2pOSzD/ARYCZwM7BuWQGamVn3K5KMPg78KiIOA+7Mtv09Ii4HtgdGAZPLCc/MzHpBkWT0bmBG9nVlsusogIh4DTgf+GSRYCS9S9JJkq6TNF9SSNq6TrtHs321j2PqtF1e0jmSZkt6XtJMSZsXic/MzNqjyACG56qOew5YCKxUtb8fWKFgPGsAnwJuBa4EdmrQ9lrg4Jptf69+IWlMdp63AV8F5gBfB66UtHlE3FYwTjMzK1GRZPQQMAEgIhZIuot06e4sSQJ2Ax4vGM+1EbEcgKRdaJyMno2IGwY5377AesDGEXFrdt4ZwD3AD0j3uczMrMOKXKa7AthdUqXcz6nAhyU9BDxAum90ZpFgImJhkeMa2BW4o5KIsvd4mXQpcZKkt5f8fmZmVkCRntExwHlkw7kj4uTscthepHtIpwPHlRbhwLaVNJ90v+o+4GTgpxFRXStvfeDqOsf+lVQ7bx3gpuodkuYO8r59hSM2M7O6ipQDmk/65V+97QTghLKCyuG3wF+Ah4FlSInwZNLlw4Oq2i0DPFPn+Geq9puZWYeNyDWJIuIrNZsulPQz4EBJP4qIx6qbNzpVnXOPa/TeWc/JvSMzsxJ1U226c0jfz/uqts2hfu9nfPZcr9dkZmZDrJuSUeV7qR4EcRfpvlGtDUj3t+5td1BmZja4bkpGe5MS0c1V2y4ENpC0YWWDpFGkuUxXRMS8oQ3RzMzqGXb3jCR9PPty0+x5K0nvAJ6PiD9I+hSwM/A74G+kS257AbsAx0fErKrTnQl8GZgm6RDSZbmvkSbpfqLt34yZmeUyaDKSdBZwakTcmL3eErgnIp5uU0y/qnl9RPb8GLAK8AjwDtLw8WWAl4E7gMm1C/9FxEuStgWOB04BxpCqO0yKiFvaFL+ZmTUpT89oMmmi643Z66uBzwI/b0dAEaFB9t9Amlib93xPkuI1M7NhKs89o9nA8lWvGyYLMzOzZuXpGf0Z+K6klYFns227SVqjwTFe6dXMzHLLk4y+TprDcyCLVnHdLXsMxCu9mplZboMmo4h4lDSibRRpaYhHSQnq4rZGZmZmPSP30O6IeAWYJekc4MaakjtmZmaFFSmUOqUdgZiZWe8qVIFB0lhJ35P012x58PnZ10dIGlt2kGZm1t2a7hlJGg/MJK0FNBuoLN09ATgM2EPSxIhwEVIzM8ulSM/oSGBt4CvAihExMSImkkrsfBlYi0VVE8zMzAZVJBntBJwRESdHxILKxohYEBGnAGeR6sSZmZnlUiQZLc+iS3P13MobKzaYmZk1VCQZPQVs1GD/RlkbMzOzXIoko0uB/SR9UdLrx0t6k6QvAPsCl5QVoJmZdb8i6xkdBkwCTga+J+m+bPtawLLAg8Dh5YRnZma9oOmeUUTMATYBjgHmkBbB25Q0zPtoYNOsjZmZWS6FVnrNlus+NHuYmZm1pFAFBjMzszI5GZmZWcc5GZmZWcc5GZmZWcc5GZmZWcc5GZmZWcc1nYwk3S/pYEkrtCMgMzPrPUV6Rq+SJrfOknSRpB2qywKZmZk1q0gFhvWAzYFzgG2Ai4HHJR0lafWS4zMzsx5QqEcTETdExOeBFYHPAY8AhwD3S7pK0qcljS4xTjMz62ItXV6LiBci4uyI2IK0+usFwNbAecATkk6UtHLrYZqZWTdr+V6PpCUk7QqcAOwJBHA1cAPwVeAeSTu3+j5mZta9CicjSWtLOh74O/AbUiXvHwITImL7iPgYqbd0H3BcGcGamVl3arpqt6R9gf2AzbJNVwCnARdHxGvVbSPiQUk/Bs5oNVAzM+teRZaQOAN4krSe0ekR8egg7e8m3UMyMzOrq0gy2h24JCIW5GkcETcBNxV4HzMz6xFF7hntSLo/VJek90k6q3hIZmbWa4oko8lAo8mtqwL7FIrGzMx6UjvK+IwllQwyMzPLJdc9o2zi6ipVm9aWtGWdpuOBA4AHWw/NzMx6Rd4BDFOAw0kTWgM4NHvUErAwa29mZpZL3mR0EfAoKdmcRZpXdH1NmwDmAzdHxONlBWhmZt0vVzKKiNuB2wEkvQf4TUTc2c7AzMysdzQ9zygivteOQMzMrHcNmowk7Z19eV5ERNXrhiLi3JYiMzOznpGnZzSVdD/oAuCVqtdqcEwATkZmZpZLnmS0DUBEvFL92szMrCyDJqOImNHotZmZWauKFEpdjKQ3AzuTJr1eGhFPlnFeMzPrDU2XA5J0nKSbq16LtKbRL4FTgTskNapdZ2Zm9gZFatN9GJhZ9XpHYEvgeODT2bZvtxiXmZn1kCLJ6N3AA1WvdwQeiYhvR8QFwE+B7YoEI+ldkk6SdJ2k+ZJC0tYDtP20pNslvSTpb5KOkTSmTrvlJZ0jabak5yXNlLR5kfjMzKw9iiSjUUD1wnrbkC7TVTwMrFgwnjWAT5HKCl05UCNJewE/A/4EfAT4AfBl0rDz6nZjsvNsBXwV2BV4DrhS0kYFYzQzs5IVGcDwOLAZcJqk9YDVgMOq9i9HSiZFXBsRywFI2gXYqbaBpCVIlwQviYgvZZuvlvRqFtOJEXFjtn1fYD1g44i4NTt+BnAPKYF9pGCcI9KJ0+/vyPseNGlCR97XzEaOIj2jC4B9JP0W+C0wD/h91f6NgIeKBBMRC3M02wxYATinZvvPSOso7V61bVfgjkoiyt7jZeB8YJKktxeJ08zMylUkGR1Nuhz2AVKlhb0jYi6ApD5Sb2bAS2wlWD97fkOh1oh4gZQE169pW6+g61+BJYB12hGgmZk1p0ih1JeB/bJHredI94teaDGuRpbJnp+ps++Zqv2VtgO1o6YtAJLmDvL+fYMFaGZmzSll0mtFdpmtv8xzNnq7nNsHajfYPjMzGyKlJqMhMid7Xqbq64rxwCM1bRfr/WTtoE6vKSLGNXrzrOfk3pGZWYkKJSNJHwC+AqxJ+mVfW8E7IqJdVRjuyp7XB14fHibprcDqwKU1bavvIVVsQBqefm+bYjQzsyYUKQe0N3AdadTaGGAW8FjNY1aJMda6AXgS+GzN9k8BbwGmVW27ENhA0oaVDZJGZW2viIh5bYzTzMxyKtIzOhS4D9g+Ip4oOR4kfTz7ctPseStJ7wCej4g/RMRrkr4NTJX0E+DXpFFxxwK/jogbqk53Jmky7DRJh5Auy30NWAn4RNmxm5lZMUWS0XuAf29HIsr8qub1EdnzY8AqABFxjqQFwMHA54HZpDJEh1cfGBEvSdqWNEn2FFJP7lZgUkTc0qb4zcysSUWS0d+A0WUHUhERjVaQrW73P8D/5GhX75KemZkNI0Umvf4U+ExWlsfMzKxlRXpGt5AGL9wk6f+ThlIvqG0UEde2GJuZmfWIIsmoutTPGSw+cVTZNveczMwslyLJaErpUZiZWU8rUpuutlq2mZlZS4oMYHidpNGS3plNJDUzMyukUDKS9F5JV5GqdM8Ctsi2LyfpSknblxijmZl1uSLlgDYEZpLqwJ1bvS8i/gksCexTSnRmZtYTivSMjgSeIC3n/W0WL5J6JfC+FuMyM7MeUiQZTQROj4j51F8PaBap9puZmVkuRZLRGBovoLdUwVjMzKxHFUlGDwEbN9i/LXB3sXDMzKwXFUlGPwc+WzNiLgAkfRP4MHBeCbGZmVmPKFKB4YfAJOBy0kqpAZwoaVlgBWA6cHJpEZqZWddrumcUEa+QktG/AS8CLwETSGsKfQvYISIWlhmkmZl1tyI9IyLiNeDE7GFmZtaSlsoBmZmZlWHQnlFW9qdZERHbFTjOzMx6UJ7LdKux+OTWscA7sq/nkqow9GWvZwPzS4nOzMx6wqCX6SJilYhYtfIAtiMNXDgJWCkixkfE0qSqCz8GXsjamJmZ5VLkntGJwJ8j4qCIeLKyMSKejIivAzfggQ1mZtaEIsloa2BGg/3XANsUCcbMzHpTkWQUwDoN9q9H/QKqZmZmdRVJRn8EDpC0t6TXl49Qsg/wxayNmZlZLkUmvX4D2BQ4GzhG0gOkntAEYHng8ayNmZlZLkXKAf0N2BA4FniWtJDe+7OvjwU2zNqYmZnlUrQcUD/wnexhZmbWEpcDMjOzjnMyMjOzjnMyMjOzjnMyMjOzjnMyMjOzjnMyMjOzjnMyMjOzjms6GUm6X9LBklZoR0BmZtZ7ivSMXgWOBmZJukjSDpLcwzIzs8KKlANaD9gcOIe0VMTFwOOSjpK0esnxmZlZDyjUo4mIGyLi88CKwOeAR4BDgPslXSXp05JGlxinmZl1sZYur0XECxFxdkRsAawNXEBafO884AlJJ0paufUwzcysm7V8r0fSEpJ2BU4A9iQtJ3E1afnxrwL3SNq51fcxM7PuVTgZSVpb0vHA34HfAJsAPwQmRMT2EfExUm/pPuC4MoI1M7Pu1PQSEpL2BfYDNss2XQGcBlwcEa9Vt42IByX9GDij1UDNzKx7FVnP6AzgSeAY4PSIeHSQ9neT7iGZmZnVVSQZ7Q5cEhEL8jSOiJuAmwq8j5mZ9Ygi94x2JN0fqkvS+ySdVTwkMzPrNUWS0WSg0eTWVYF9CkVjZmY9qR1lfMaSSgaZmZnlkuueUTZxdZWqTWtL2rJO0/HAAcCDrYdmZma9Iu8AhinA4aQJrQEcmj1qCViYtW8bSVuTJtbWs05E3FvVdhLwfeBfgeeAC4GDI2JuO2M0M7P88iaji4BHScnmLNK8outr2gQwH7g5Ih4vK8BBHAxcW7Pt0coXWdL6PSn+7wIrAccC60uaGBELhyZMMzNrJFcyiojbgdsBJL0H+E1E3NnOwHK6PyJuaLD/OOBOYM9K4pH0D+CPwB7AL9ofopmZDabIEhLfGyaJqCFJ7wQ2Bc6r7gFFxHRSCaPdOxWbmZm90aA9I0l7Z1+eFxFR9bqhiDi3pcjyOVXSr4HngZnA4RFxS7Zv/ey5XuK8o2r/G0ga7F5SX5FAzcxsYHku000l3Q+6AHil6rUaHBNAO5NRP/Aj4BrgGWAd4NvAnyRtFRE3AstkbZ+pc/wzwHvbGJ+ZmTUhTzLaBiAiXql+3UkRcRtwW9WmmZIuIfWCjgK2r24+0GkGOPe4Ru+d9ZzcOzIzK9GgySgiZjR6PVxExJOS/gjslG2akz0vU6f5eOr3mMzMrAPaUYGhk97Eoh7PXdlzvXtDG1D/XpKZmXVAMwMYmjJEAxheJ2kFYBJphVki4m+S/gJ8RtKPqoZ2bwe8E5g2lPGZmdnAmhnA0GjAQq22DmCQ9DPgYeBW4FnSirIHA0sCh1Q1PZg0p+h8SaexaNLrjcCv2hWfmZk1J/cAhmHmDuCTwFdJhVnnkEbW/Wf1HKiIuErSDsD3gN+RygFdBHwr73pMZmbWfk0PYBgOIuIY0kqzedpeBlzW3ojMzKwVRVZ6NWvKidPv79h7HzRpQsfe28zyG+kVGMzMrAuM1AoMZmbWRUZkBQYzM+suXVOBwczMRq7CAxgkjQa2BlbLNj0MzIiIl0qIy8zMekihZJQNYjgBWJpF944CmCvpmxExtZzwzMysFzSdjCTtSRrEMAv4IXA3KSGtC+wPnCnpxYjwKqpmZpZLkZ7Rd4B7gc0iYl7V9oslnUwqtXMoXtLbzMxyKlK1ey3g7JpEBEBE9ANnA2u2GpiZmfWOIsnoSRrPMVoIPFUsHDMz60VFktFUYLKkt9XukLQUsC+pd2RmZpZLnnJAW9ZsuhbYAbgju0d0L2kk3brAAcBsYGbJcZqZWRfLM4DhGhatnlpRuUx3bNW+yrb3ANOBJVoNzszMekOeZDSl7VGYmVlPy1MO6JyhCMTMzHpXkQEMZmZmpWqlNt0mwPtJJYFqk1pExPdbCczMzHpHkXJASwLTgA+SBi1Ur20UVducjMzMLJcil+kOIyWio0hrGwnYB/gIaUj3zaRh3mZmZrkUSUYfB34VEYcBd2bb/h4RlwPbA6OAyeWEZ2ZmvaBIMno3UFlgb0H2PAogIl4Dzgc+2XpoZmbWK4oko+dYdK/pOVItupWq9vcDK7QYl5mZ9ZAiyeghYAJARCwA7iJdukOSgN2Ax8sK0MzMul+RZHQFsLukSrmfU4EPS3oIeIB03+jMkuIzM7MeUGSe0THAeWTDuSPiZEljgL1I95BOB44rLUIzM+t6iqitgWqNSJrb19fXN3fu3ELHnzj9/pIjsuHooEkTOh2C2bAybtw4+vv7+yNiXL39LgdkZmYd10o5oNHA1sBq2aaHgRkR8VIJcZmZWQ8plIwk7Q2cQKpLV10KaK6kb0bE1HLCMzOzXlCkNt2epKXHZwE/BO4mJaR1gf2BMyW9GBG/KDFOMzPrYkV6Rt8hLTW+WUTMq9p+cbYM+Y3AoYCTkZmZ5VJkAMNawNk1iQiAiOgHzgbWbDUwMzPrHUWS0ZMsuk9Uz0LgqWLhmJlZLyqSjKYCkyW9rXaHpKWAfUm9IzMzs1wGvWckacuaTdcCOwB3ZPeI7iWNpFsXOACYTVrXyMzMLJc8AxiuISWbapXLdMdW7atsew8wHVgCMzOzHPIkoyltj8LMzHraoMkoIs4ZikDMzKx3uTadmZl1XNFyQGOBbwG78sbadNOA4yPi+XLCMxuZOlWd3dXCbaQqUg5oPGm03DqkkXO3ZbsmAIcBe0iaGBHPlBalmZl1tSKX6Y4E1ga+AqwYERMjYiKwEvBlUoWGI0qL0MzMul6RZLQTcEZEnBwRCyobI2JBRJwCnAXsUlaAZmbW/Yoko+VZdGmunluzNmZmZrkUSUZPARs12L8Rrk1nZmZNKJKMLgX2k/RFSa8fL+lNkr5Aqk13SVkBmplZ91NEbaWfQQ6QlgGuB1YHngbuy3atBSwLPAhsHhFzSoxz2JA0t6+vr2/u3LmFju/UkF+zdvOwcmtk3Lhx9Pf390fEuHr7m+4ZZUlmE+AYYA6wafaYDRwNbDqcEpGkt0n6saR/SHpR0l8k7dTpuMzMbJGm5hlJWhLYA7gvIg4lreg63F0IvJc0SfcRYDJwoaQdI+L3nQzMzFrXyasN7g2Wp9lJry8DZwAHkpYXH9YkfRTYHtgtIi7Mtl1NqhrxX4CTkVlJfAnaWtHUZbqIWAjMApZqTzil2xXoBy6ubIh0k+wcYG1J63YqMDMzW6TIAIb/AD4BbBIRL7clqpJIup6Ufzav2f5+4AZgz4j4Zc2+wUYm9AH09fUViunl1xYWOs7MbDgY/eZi9bX7+/sh/T6ue4IihVL/DOwG/G+20usDwAu1jSLi2gLnLtsyQL1rB89U7S8i+vv75xU4rpLB+gu+r3WGP7eRyZ9bG7xU/NClgAH/Gi+SjKZXfX0S9VeBDYbPSq+Nun6L7Rto2GEZKr2udr6Hlc+f28jkz21kKZKM9qXxL/jhZA71ez/js2dXFjczGwaaTkYRMbUNcbTLXcDukt6UDb6o2CB7vrMDMZmZWY2m7kRJWlbS+yWt3q6ASnYhMA7YsWb73qS5UncPfUhmZlYrV88oq0F3MvA50j2hyki1XSPi6faF17LfA1cDZ2ZljB4B9gG2AHbuZGBmZrZI3p7RV4AvAE+Slha/A9gcOLVNcZUim1O0C3AB8APgD8C/kCbBXtrJ2MzMbJFc84wk/QVYEtgsIp7Ltp1OKq2zbEQUqxraYzy6Z2Ty5zYy+XMbWfL2jNYCplYSUea/ScO3XZzJzMxaknc03VjgiZptT1Ttsxz8F9rI5M9tZPLnNrI0M5qu9npe5bVKisXMzHpUM/OMPippharXbyUlpD0kbVjTNiLixJajMzOznpB3AEOz1T0jIoZLOSAzMxvm8l6m26bJx7alRzpMtLJyrKTVJV0kqV/Sc5J+P9AyFpIOlHS/pJclPSTpW9l8LytgKD43SUdJ+oOkf0oKSUeU/o30mHZ/bpImSDpB0m1ZuzmSZno16A6ICD+aeJAKxc4B9iMl3XOBBcBHBzluOdKgj/8lzX3agbSMxdPAu2rafjc755HA1sB3gFeBYzr9/Y/UxxB9bvOBPwGnkS5hH9Hp73ukP9r9uZHmUN6T/R+bBHwUOD/7/L7e6e+/lx4dD2AkPbJ/qEGqPFHZJuA64J5Bjj0OeBFYqWrbMsA84JSabS8CJ9Ucf1SWkN7V6vfRa4+h+Nyy7W/Knsc5GY2Mzw14B9ntiprjrwZmd/pn0EuPQS/7SNpusDYNjt2+6LHDVCsrx+4KTI+I14fIR8Qc4FLS+lAVHwbGZOesNpU04MSXD5o3FJ8b8cZivNa6tn9uETE7O2etm4FlJC3Z2rdgeeW5B3GZpKsk7SBp0EEJkt4iaVdJM0i14brJ+sDddX7p/LVq/2Kyf9CrU79K+F+B5SQtV3WOIFUcf11EPED6S6/ue1hDQ/G5Wfk68rlJEune98MR8QUEbmkAAAwzSURBVGLTUVsheYZ2bwScAFwCzJY0HbgJeIi0HpBI6wOtCWxGuq67NPBHoHbI90hXdOXYpUk/p3rrJ1Uf+8/s+YWov6T7sw3ewwY2FJ+bla9Tn9vXgE1Ia7fZEBk0GUXEncAHJX0A+BKp2vWnqL/C6zxSIdVTIuLmkmMdLppaObbgsa28h9U3FJ+blW9IPzdJuwA/JJU/O3uQ81uJck96jYjrgeuzS3UbA+sCy5I+1KdJXeLbuvy6edGVY58l/ZzyHDsHGCtpdJ3e0dIN3sMGNhSfm5VvSD83SR8DfkH6g/pzTUVqLSuy0usC0mW6m8oPZ9grtHJsRLwo6WHqX+PeAHg6IiqXDO4i9TLXA26tNJK0Bqlyulenbd5QfG5WviH73CR9hJSE/gB8Jvs9Z0PIkyib08rKsRcCk6pLKkkan51rWlW7PwAvA5+tOX4f4DXSaCBrzlB8bla+IfncJH0oa38F8ImIeLWE2K1ZnR5bPpIepB7LVcBs0s3NbUhDrhcCO1a1u4ZsFGrVtuVJixPeSrrv9jHgetKliJVr2h5OSjxHAFsB3wZeAY7v9M9gJD6G8HPbCvg46ZdlAL/MXn8ceGunfw4j7TEUnxtp1ecXSAOytiINwqp+jO70z6FXHh0PYKQ9gKWAn2T/0F/K/rHvUtNmsf8c2fY1SXMm5pFm6/8BWK9OOwFfBx4g9ZIeAQ4hm1Tpx7D93K7JklC9xyqd/hmMxEe7PzfSH3wDfWb+3IbwkatQqpmZWTv5npGZmXWck5GZmXWck5GZmXVc0/OMIK0BQpoHsxxVk14j1U8zMzNrSu4BDJLWAfYH9iANm4Q06gsWldZ4ijSc9dSIuKfEOM3MrIsNmowkrQ4cSyrJ/iIwkzRe/yHSmP1KodQ1SOPyJ5IqBUwDDo6Ih9sVvJmZdYc8yehl4A7gJGBaRDw/SPuxpEl+B5LG9I8pKVYzM+tSeZLRzhFxccNGbTjWzMx6x6Cj6SLiYklbFDm5E5F1G0nrS3pN0qROx9IsSVtLCkmTG20bolgelXRNC8fvIukVSWuWGJZ1UN6h3dMl7dnWSMzaRNLHsl+4R9TZ9zZJt0p6WdKWOU53AvCniJg+wHstLeml7P32ajH0jpC0oaQjJK3S6VgGEhEXkW4fHNvpWKwceZPRPcDPJH0rT2NJ2xcPyaxcEfE74H+BAyW9vbI9W5vrAtKKxFMi4tpG58kWmJxESkgD+QwwilRPcL8WQx8K15IGHJ1XtW1DUrHeVToRUBNOAnaVtF6nA7HW5U1GE4HLgKMlnSKp7nGSdpJ0I3B5WQGaleQo0uKEB1Rt+29SNedDI+LnOc7xJdII0t83aLMfcDXwI2CrbDTqsBURCyPipRiZ6/dMI1Xc3r/TgVjrciWjbATdTsCpwBeBS7JRcyj5pKTbSWuCrAn8Z5viNStqGnAv8A1JYyT9OykxnR4RRw92sKQ3A7sA02OA9W4kvZfUqzgH+BnwKjBlgLaTs0t520k6TNJjkl6UdKOkzbI2W0m6TtLzkv4h6T8GOMf22WW1x7LLjX+V9Mk8P5Tae0bZpczKcttXZ/tC0tTK/uz1KnXOtdh9IEnvlvRLSf2S5km6tFGCljRa0nck3ZVd7pybHbNRbduImE+aarJHnu/Vhrdmlh1fCHxJ0iPAMcAMSWcCBwGrk/5i/C7wk4h4rh3BmhUVEQslHU1KFBeQ/ri6jNTbyWNj4G00XuF4P+B54DcR8byk3wH7SDos3rhSabVjgCVIl5xGAd8ELpe0D3AmcBopsX0COFLSIxHxPzXnOBYYC5xCmoA+BThf0piImJrz+6uYBqwIfAH4AekSPaR5hU2RNI50GfDdwE+Bu0lrBl1NujRY2/4tpM9kc9Jlw58AfcDngT9J2jIi/lJz2PXAhyStHRH3NhujDSPNrjkBjCb1gBYCC4AnSAlpyU6vh+GHH40epD++HiH9wr4NeHsTx07JjttpgP1jgGeAqVXbds6O+Uid9pOzfbcCo6q275Rtfw3YtGr7KOAfwPV1zvEY0Fe1vS/b9kz1/0tg66z95EG2Vc67dZ24j2CAdX6AR4Frql7/IGs7pabdj7Lt19RsPyjb/qGa7UsBs2rbZ/v2yo7ZvdP/vvxo7ZG7UKqksZL+LfvPvFP2D29O9p/kpoh4Me+5zDpkVVLvBuCMaK4Hv2z2/MwA+3cj3ZM6p2rb74B/klYpHcgpEfFK1euZ2fMNEXFzZWPW5ibSZfB65+ivattP6oksTUo2nbILqUTYuTXbBxoBtxfpUuotkt5ReZB+x0wHtpBU26Oakz0vV1LM1iG5kpGkw0h/aR0HzCX95bQm8P+y11d46LcNZ5KWJa30uQTpF+RB2Wi6vCqzwzXA/v1IBYP/JmkNSWuQRqNNB3bKfqnW84ZyWRHxbPblI3XaPgssU2d7vTqQd2fPqw3wvkNhNeCBqBkcERH/IP3eqLUOsDbp51j72Jf02dX+HGvrY9oIlfee0RGkobFfJJUEqnzwD2Q3W39HGvr9nog4rvwwzYqT9Fbgt6R7F5OATYD/AvYE8oyig/QLEVIdxtrzrwpsQ/rFeP8Ax+9FujxVa6BRbM2Mbqv3i3igpNmqRr/06/0+Gah9vfhEmjv0jQbv8XTN6/EDbLcRJm8y2iEi6g5njYjZkrYm3RQ+OvuP+eUY+Iat2ZDJpiGcD2wK7BUR10r6C3AIcIik86v+uGrkzuy53mWyKaRfpJ+n/l/8/0nqOdVLRmVYF7ikZts62XORQsWNfh6Vy5TjSZfqAZA0hjTw4cGqtg8DEyQtUd07krQi6b5WrQdIl0OvauL3xxrZ850NW9mwl3dod6N5FWT3i3YhXaf+Iov/xzDrlP8m3eN8fS5RRLwAnAisTxpkkMdtwDxSZfrXZcluMnBHRJwREb+ufZCS4fqSNi3lO1rcAZJe/+Wefb0/KTHOKHC++dnzYr1AFvX8aie2H8Tiv08uJi03s3fN9oMHeN9zgRUYoGckafk6mzcDnoqI+wY4p40QhRbXqyf76/LLkh4jjaIx66isYsiXqD+X6CfAvwOHAhcNdq6IWCBpGrCzpNER8XK264Oky39nNjj8N6RL3fsBNzdoV9Rs4EZJZ5F6aFOAlYHPZYm3WTeTRsseKmlp0nD1RyLiRuAK0iCDIyUtQ7q3tQUpKcyuOc9xwKeB0yVtDNxFGlDxgTptIQ1vnwQcL2lb4CrSHwArA9sBL5EuhwKplBNpQv5ZBb5HG2ZKX3Y8u2f0mbLPa9aMbNLnMQwwlygi5pF6TZtI+mDO055CGqG2Q9W2SsmfaQMdFBF3knoUn6wzGqwMBwO/AL4CHEkaFv6ZiGiUIAcUEbNIAwaWJH3P55NVrsgut+0MXAN8lfQzHkWaP/R8zXmeJSWLi0i9o+OAt5ISymJL0USaTPwx4Guky3XfI/Vg9yRd8qv9g2L37HynFvk+bXjJs4TEdhFxZaGTS9tHxBWFIjMbhiRdBoyNiInDIJbJpGoJ20TENZ2NZuhJugV4LCJ263Qs1ro8PaPLJF0laYc8Q2ElvUXSrpJm0LiGl9lI9E3gA030pqwNJO0CbMDA959shMlzz2gjUpXiS4DZkqaTJt89RBpZU1l2fE3SdeNts9eXk+p0mXWNiLiLEu+1WjGRlpAY1ek4rDyD/qfKrnd/UKl8/pdI14s/xeLDP0W62TiNNCO8HTdqzcysCw16z2ixA9Kluo1JcxuWJSWlp0nj/G/z/CIzM2tW08nIzMysbIWufWeTzz5Bqr01n1R5eHrBOQ1mZtbjilymm0gaJfdW3lhfag7w/Yj4cXnhmZlZLygy6fWH2fO+pJnRqwCfJI2u+5GkvIUnzczMgGI9o+eB/4qIw+rs24+0MuXXIuIn5YRoZmbdrkjP6DnSqouLycqPXEAq0mhmZpZLkWR0NfDRQfavXiwcMzPrRUWS0WnA5pK+NsD+VYAnCkdkZmY9p8g9o4WkqsBLADOBM4Bbsm1bkarsHhwRJ5cbqpmZdasi84y+D/wrqe7cltmjOqPdBsyVtD5wb0S81nKUZmbW1VqqwCBpHKmQ6obZYyNgbVKSC+BV4D7grxHx2ZajNTOzrlR6OSBJo0jLOVeS00bA+hExrtQ3MjOzruHadGZm1nGlLztuZmbWLCcjMzPrOCcjMzPrOCcjMzPrOCcjMzPruP8D0joy6/UecG4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the angular velocity\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(omegas, bins=10, alpha=0.5, density=True)\n",
    "ax.set_xlabel('$\\omega$ (Angular velocity)')\n",
    "ax.set_ylabel('$p(\\omega)$ (Probability density of ang. vel.)')\n",
    "\n",
    "# Plot the amplitude\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(Xs, bins=10, alpha=0.5, density=True)\n",
    "ax.set_xlabel('$X$ (Amplitude)')\n",
    "ax.set_ylabel('$p(X)$ (Probability density of amplitude)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "+ What does the probability density in the figures above represent? Is the uncertainty aleatory or epistemic?\n",
    "+ Rerun the code above giving different values to ``num_samples``. Can you trust the results when you pick small values? How can you pick the right value for ``num_samples``?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The model calibration problem\n",
    "\n",
    "The model calibration problem is the inverse of the uncertainty propagation problem.\n",
    "That is why such problems are also called **inverse problems**.\n",
    "It goes as follows. \n",
    "One observes a quantity that is predicted by the model and they want to go back and characterize how this observation changes the state of knowledge about the parameters of the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Driving a trailer on a rough road (2/3)\n",
    "\n",
    "In this example, imagine that we put sensors on the suspension of the trailer to measure the amplitude of oscillation.\n",
    "Now, the amplitude we measure is not exactly the prediction of the model.\n",
    "Why is that?\n",
    "First, there is observation noise.\n",
    "Second, our model has a systematic error (e.g., we have completely ignored any dumping effects).\n",
    "So, *the measurement is not the same as the model prediction*.\n",
    "This means that we need to add one more node to the graphical representation of the model.\n",
    "The new node, let's call it $X_m$ is the result of the measurement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: Trailer Pages: 1 -->\n",
       "<svg width=\"314pt\" height=\"260pt\"\n",
       " viewBox=\"0.00 0.00 314.00 260.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 256)\">\n",
       "<title>Trailer</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-256 310,-256 310,4 -4,4\"/>\n",
       "<!-- k -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>k</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\">k</text>\n",
       "</g>\n",
       "<!-- X -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>X</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"135\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"135\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">X</text>\n",
       "</g>\n",
       "<!-- k&#45;&gt;X -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>k&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M45.81,-148.81C63,-137.67 88.62,-121.06 107.99,-108.5\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"109.92,-111.43 116.4,-103.05 106.11,-105.56 109.92,-111.43\"/>\n",
       "</g>\n",
       "<!-- m -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>m</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\">m</text>\n",
       "</g>\n",
       "<!-- m&#45;&gt;X -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>m&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M107.35,-144.76C111.71,-136.28 117.15,-125.71 122.04,-116.2\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"125.23,-117.64 126.7,-107.15 119.01,-114.44 125.23,-117.64\"/>\n",
       "</g>\n",
       "<!-- y0 -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>y0</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"171\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"164\" y=\"-158.8\" font-family=\"Times,serif\" font-size=\"14.00\">y</text>\n",
       "<text text-anchor=\"start\" x=\"171\" y=\"-158.8\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\">0</text>\n",
       "</g>\n",
       "<!-- y0&#45;&gt;X -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>y0&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M162.65,-144.76C158.29,-136.28 152.85,-125.71 147.96,-116.2\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"150.99,-114.44 143.3,-107.15 144.77,-117.64 150.99,-114.44\"/>\n",
       "</g>\n",
       "<!-- omega -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>omega</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"243\" cy=\"-162\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"237.94\" y=\"-157.8\" font-family=\"Times,serif\" font-size=\"14.00\">ω</text>\n",
       "</g>\n",
       "<!-- omega&#45;&gt;X -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>omega&#45;&gt;X</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M224.19,-148.81C207,-137.67 181.38,-121.06 162.01,-108.5\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"163.89,-105.56 153.6,-103.05 160.08,-111.43 163.89,-105.56\"/>\n",
       "</g>\n",
       "<!-- v -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>v</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"207\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"207\" y=\"-229.8\" font-family=\"Times,serif\" font-size=\"14.00\">v</text>\n",
       "</g>\n",
       "<!-- v&#45;&gt;omega -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>v&#45;&gt;omega</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M215.35,-216.76C219.71,-208.28 225.15,-197.71 230.04,-188.2\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"233.23,-189.64 234.7,-179.15 227.01,-186.44 233.23,-189.64\"/>\n",
       "</g>\n",
       "<!-- L -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>L</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"279\" cy=\"-234\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"279\" y=\"-229.8\" font-family=\"Times,serif\" font-size=\"14.00\">L</text>\n",
       "</g>\n",
       "<!-- L&#45;&gt;omega -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>L&#45;&gt;omega</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M270.65,-216.76C266.29,-208.28 260.85,-197.71 255.96,-188.2\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"258.99,-186.44 251.3,-179.15 252.77,-189.64 258.99,-186.44\"/>\n",
       "</g>\n",
       "<!-- Xm -->\n",
       "<g id=\"node8\" class=\"node\">\n",
       "<title>Xm</title>\n",
       "<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"135\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"124.5\" y=\"-14.8\" font-family=\"Times,serif\" font-size=\"14.00\">X</text>\n",
       "<text text-anchor=\"start\" x=\"134.61\" y=\"-14.8\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\">m</text>\n",
       "</g>\n",
       "<!-- X&#45;&gt;Xm -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>X&#45;&gt;Xm</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M135,-71.7C135,-63.98 135,-54.71 135,-46.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"138.5,-46.1 135,-36.1 131.5,-46.1 138.5,-46.1\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a2642aed0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g = Digraph('Trailer')\n",
    "g.node('k')\n",
    "g.node('m')\n",
    "g.node('y0', label='<y<sub>0</sub>>')\n",
    "g.node('omega', label='<&omega;>')\n",
    "g.node('v')\n",
    "g.node('L')\n",
    "g.node('X')\n",
    "g.edge('v', 'omega')\n",
    "g.edge('L', 'omega')\n",
    "g.edge('y0', 'X')\n",
    "g.edge('omega', 'X')\n",
    "g.edge('k', 'X')\n",
    "g.edge('m', 'X')\n",
    "g.node('Xm', label='<X<sub>m</sub>>', style='filled')\n",
    "g.edge('X', 'Xm')\n",
    "g.render('trailer_m_g', format='png')\n",
    "g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have filled the node $X_m$ with color to indicate that it is observed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving inverse problems\n",
    "\n",
    "We will need a couple of lectures to understand what is the right way to pose and solve the problem.\n",
    "But here is the answer:\n",
    "+ Quantify our **prior** state of knowledge about all the model parameters (by assigning probability densities to them).\n",
    "+ Use Bayes' rule to condition the prior knowledge on the observations. This updated knowledge is our **posterior knowledge**. Unfortunately, this posterior knowledge is rarely analytically available. This is why we need the third step.\n",
    "+ Create a practical procedure that characterizes our posterior state of knowledge.\n",
    "\n",
    "The majority of the lectures of this class are about the third step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hands-on activity: Catalytic Conversion of Nitrate to Nitrogen\n",
    "\n",
    "The purpose of this hands-on activity is to:\n",
    "\n",
    "+ Familiarize yourself more with the Python language and Jupyter notebooks.\n",
    "+ Experience a real model calibration problem and appreciate its difficulties.\n",
    "+ Gain some additional experince with the uncertainty propagation problem using Monte Carlo sampling.\n",
    "\n",
    "This is Example 3.1 of [(Tsilifis, 2014)](http://arxiv.org/abs/1410.5522).\n",
    "\n",
    "Consider the catalytic\n",
    "conversion of nitrate ($\\mbox{NO}_3^-$) to nitrogen ($\\mbox{N}_2$) and other\n",
    "by-products by electrochemical means.\n",
    "The mechanism that is followed is complex and not well understood.\n",
    "The experiment of [(Katsounaros, 2012)](http://www.sciencedirect.com/science/article/pii/S0013468612005208) confirmed the\n",
    "production of nitrogen ($\\mbox{N}_2$), ammonia\n",
    "($\\mbox{NH}_3$), and nitrous oxide ($\\mbox{N}_2\\mbox{O}$) as final products\n",
    "of the reaction, as well as the intermediate production of nitrite ($\\mbox{NO}_2^-$).\n",
    "The data are reproduced in [Comma-separated values](https://en.wikipedia.org/wiki/Comma-separated_values) (CSV) and stored in\n",
    "[data/catalysis.csv](data/catalysis.csv).\n",
    "The time is measured in minutes and the conentrations are measured in $\\mbox{mmol}\\cdot\\mbox{L}^{-1}$.\n",
    "Let's load the data into this notebook using the [Pandas](http://pandas.pydata.org) Python module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Time</th>\n",
       "      <th>NO3</th>\n",
       "      <th>NO2</th>\n",
       "      <th>N2</th>\n",
       "      <th>NH3</th>\n",
       "      <th>N2O</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>500.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>30</td>\n",
       "      <td>250.95</td>\n",
       "      <td>107.32</td>\n",
       "      <td>18.51</td>\n",
       "      <td>3.33</td>\n",
       "      <td>4.98</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>60</td>\n",
       "      <td>123.66</td>\n",
       "      <td>132.33</td>\n",
       "      <td>74.85</td>\n",
       "      <td>7.34</td>\n",
       "      <td>20.14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>90</td>\n",
       "      <td>84.47</td>\n",
       "      <td>98.81</td>\n",
       "      <td>166.19</td>\n",
       "      <td>13.14</td>\n",
       "      <td>42.10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>120</td>\n",
       "      <td>30.24</td>\n",
       "      <td>38.74</td>\n",
       "      <td>249.78</td>\n",
       "      <td>19.54</td>\n",
       "      <td>55.98</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>150</td>\n",
       "      <td>27.94</td>\n",
       "      <td>10.42</td>\n",
       "      <td>292.32</td>\n",
       "      <td>24.07</td>\n",
       "      <td>60.65</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>180</td>\n",
       "      <td>13.54</td>\n",
       "      <td>6.11</td>\n",
       "      <td>309.50</td>\n",
       "      <td>27.26</td>\n",
       "      <td>62.54</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Time     NO3     NO2      N2    NH3    N2O\n",
       "0     0  500.00    0.00    0.00   0.00   0.00\n",
       "1    30  250.95  107.32   18.51   3.33   4.98\n",
       "2    60  123.66  132.33   74.85   7.34  20.14\n",
       "3    90   84.47   98.81  166.19  13.14  42.10\n",
       "4   120   30.24   38.74  249.78  19.54  55.98\n",
       "5   150   27.94   10.42  292.32  24.07  60.65\n",
       "6   180   13.54    6.11  309.50  27.26  62.54"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import io\n",
    "import requests\n",
    "url=\"https://raw.githubusercontent.com/PredictiveScienceLab/uq-course/master/lectures/catalysis.csv\"\n",
    "s=requests.get(url).content\n",
    "catalysis_data = pd.read_csv(io.StringIO(s.decode('utf-8')))\n",
    "catalysis_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the data using [Matplotlib](http://matplotlib.org):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x1a26d21310>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEZCAYAAABiu9n+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydeVxV1fr/3+sAMnsYRQUFxRHJMRVFQwu42oC3rK9Zmea9Dd977ec3UxEbvHbL4aZlZrc0x2uOXVMrR8qpcsB5VhRBBREFBBRBpvX7Yx+Q4aDMh2G9X6/9Opy1nr32sw+wP+dZzxqElBKFQqFQKIyhM7UDCoVCoai9KJFQKBQKRakokVAoFApFqSiRUCgUCkWpKJFQKBQKRamYm9qBqkIIkYMmemmm9kWhUCjqEI2BPCmlUT0Q9WUIrBAiDxB6vd7UrigUCkWdITU1FUBKKY32LNWbSAJI0+v1+pSUFFP7oVAoFHUGBwcHUlNTS+2BUTkJhUKhUJSKEgmFQqFQlMpDRUIIMUAIIUs5OhSzDRJC7BdCZAghbggh5gshHIy0aSeEmCuEiDfYHhJChFTljSkUCoWi8pQnJxEK7ClWFpP/gxBiALAZ2AC8DzQHZgK+Qoj+Usq8QuetB7oDE4FoYBSwXgjxjJRyc/luQaFQKBTVRXlEIlJKuf8B9f8CTgHD8gVBCBEPbAdeANYYyp4EAoHnpJTrDWU7gdbAbDShqXa+2R1FZw89fb1djNbvjUrkRGwqbwV414Q7CoVCUSupktFNQgh3oCfwbuGIQUoZLoSIA4ZiEAngWSAV2FjITgohlgELhBA+UsozVeHXg+jsoWfMyqPMe6lbCaHYG5VYUKeoPUgpSUxMJDMzk7y8vIefoKgUOp0OKysrXFxcEEKY2h2FiShP4nq+ECJHCJEqhPhZCNGjUJ2v4fWUkfNOFqrPtz1TrPsJ4ESxtqqVvt4uzHupG2NWHmXzyWv88+cz3Lx9r4hAlBZlKGoeKSVxcXEkJiaSnZ1tancaBNnZ2SQmJhIXF0d9mU+lKD9liSRSgTnALiAZ6AhMAv4QQgRIKQ8AzgbbZCPnJ6PlH/JxBiJLscuvL4EQ4mETIMo9i66vtwtzX+zGq4sPkCfh2q0MDsQkK4GohSQmJnL79m3c3NxwcnIytTsNhuTkZBISEkhMTMTV1dXU7ihMwEMjCSnlUSnlO1LKjVLK36SUC4C+wB3gk+LmpTVTRruH1VU5/dq6MKxnCwC2nL7OP4d0UgJRC8nMzMTS0lIJRA3j5OSEpaUlmZmZpnZFYSIqNE9CSnkdLSHtZyhKMrwaiwKcKBphJD3ADoxHI0gpHR50oEU8FSLYx63g5/VH4yrajKIaycvLw8zMzNRuNEjMzMxUDqgBU5nJdDruf+s/bXg1lk94hKK5itNARyFE8Ws/Yng1lteoNvZGJfLu9yd4tY8nAL+cvcFGJRQKhUIBVFAkhBBNgSBgP4CUMhY4BLxc+OEvhHgCcAd+KHT6esABeKZYs68C52tiZFM+hZPUk5/sSDO9FQAT1p1gb1RiTbmhUCgUtZayzLheIYT4pxDiWcPs67eAfYA1EFbINBToAqwSQjwhhBgBLAcOAN8XstsM7AQWCSFGCyEGCiGWAv2ACVVyV2Wg+CgmKwsz3n68LQDZOXn873dHlFAoFIoGT1kiiZNo3/qXAOHAP9Ae/L2klIfyjaSUO4CnAS9gE/CZ4XWwlDK3kJ0E/gysBqYBW4DOaJPrfqr0HZWRE7GpJUYxvfCoBy2dbJBAp+aNORFb4TSHQlFmli5dihACa2trYmNjS9R37dqVAQMGFClLTExk4sSJdOjQAWtraxwcHAgICGD58uUlhqsmJiby/PPP06ZNG+zs7GjcuDHdunVj3rx55ObmolA8iIcOgZVSzgBmlKUxKeVWYGsZ7NKAMYbDJBibSW1hpmPsE2159/vj7LuUxIfP+JjAM0VDJTMzkylTprBo0aIH2p07d47AwEAyMzOZMGECvXr1Ij09ne+//56RI0eyadMmVq5ciU6nK2jXysqKsLAwPD09ycnJYcuWLbz99tucOnWKb775piZuT1FXkVLWiwNI0ev1srLk5ObJx2ftlJ6hP8s3/nOw0u0pqoaYmBgZExNTZe19veui/OPizVLr/7h4U36962KVXe9BLFmyRAJy0KBB0szMTJ45c6ZIfZcuXWRAQICUUsqcnBzp4+MjHRwc5MWLJf2bNWuWBOT06dMfet1hw4bJRo0ayezs7AfaVfVnr6hd6PV6CaTIUp6taqnwYpjpBP8X2A6AbacTOKm6nOol+cuyGMs75eerOnvU7C6HoaGh6PV6Jk+eXKrN+vXrOXPmDGFhYXh7l4yGx40bR8eOHfn0008fOjPdxcUFnU5XEHEoFMaoTzvTVRlPPdKMr3Ze5Nz123wWfp4lr/UytUuKUsjJzSM+tfwTvVo42jDlaR/+97sjfBTSie6ejgAcuXyLD388zUchnWjhaMPV5LvlareZ3gpzs4o9dB0cHAgLC2PChAns27ePPn36lLAJDw8HICTE+Mr6QghCQkKYOXMmhw8fxs/Pr6BOSklubi63b98mPDycpUuXMmHCBCUSigeiRMIIOp1gXFA73lh+mJ3nb3L4cjI9PNVM39pIfGom/f+1s1JtjF1zrExlZeG3iQNp4WRTYV/GjBnD3LlzmTRpErt37y5Rf+XKFQBatWpVahteXl4FtoVF4quvvuLtt98GNDGZPHkyH330UYV9VTQM1FeIUgjycSvobpi93dhSUwpF1WNlZcXUqVPZs2cPmzZtqlAb0jC6qfjKrcOGDePgwYNs376dSZMmMWvWrALRUChKQ0USpSCE4N3g9oxcHMHeqCT2Xkykbxu1plNto5neit8mDqxUG0cu3yqIHL4Y1rWg66mi/lSWkSNHMnv2bMLCwhg8eHCRupYtWwIQHR1Nhw4djJ3O5cuXAWjRokWRcldX14JF+oKCgnB2dmb8+PGMHj2abt3UsvgK4yiReACPtXWhp5cjB2NuMTs8kj7ezmpd/VqGuZmuUt07e6MSmfrzGVa+3hugViwTr9PpmDZtGkOGDGHFihVF6gIDA1mwYAEbN240KhJSSn788UecnJzo3r17ifrC9Oql5doiIyOVSChKRXU3PYD8aALg8OVb7Iq8aWKPFFVJ8Vn3hfcYMfVs+5CQEPz9/fnwww/JysoqKH/22Wfp0KEDM2bMICoqqsR5n332GWfPnmX8+PE0atTogdfYuVPL5bRp06ZqnVfUK1Qk8RD8WjvTr40Lv19MZPb28wxo56qiiXpAaZtLFRYKU0cUM2fOpF+/fgA0adIEAHNzc9atW0dgYCB+fn4lJtP95z//YejQoUycOLGgnVmzZnH27Fkef/xx3N3dSU1NZfv27cyfP5+hQ4fSo0cPo9dXKEBFEmViXLA2b+JUXBrbTieY2BtFVWBsWZZ88oXC1Muy+Pv7Gx3q6uPjw7Fjxxg5ciQLFy5k0KBBDB8+nKioKBYvXszatWuLLKvetWtX4uLiGD9+PMHBwbzyyitEREQwe/ZsVq9eXZO3pKiDiPyREHUdIUSKXq/Xp6Q8bAO7ijF66UF2nLtBezd7toztj06noomaJD8Z6+npaWJPGh7qs6/fODg4kJqamiq1fXlKoCKJMjIuSIsmzifc5ueT8Sb2RqFQKGoGJRJlxNddz2DfpgDMCY8kJ1ft1KVQKOo/SiTKwTtB7RACLiWmq21OFQpFg0CJRDlo52bPkC7NAfji1wtk5ahoQqFQ1G+USJSTsYHtMNMJYm9lsPbQVVO7o1AoFNWKEoly0srFlue7ewDw5Y4LZGarnb0UCkX9RYlEBXj7iTZYmAkS0u6x4sAVU7ujUCgU1YYSiQrg4WjD8F7aQmtf77rI3awcE3ukUCgU1YMSiQry94FtsDTXkXgni2V7L5vaHYVCoagWlEhUELfGVozw02agfrM7irTMB28VqVAoFHURJRKV4K0B3tg0MiM1I5vFv0eb2h1FHWXp0qUIIbC2tiY2NrZEfdeuXRkwYECRssTERCZOnEiHDh2wtrbGwcGBgIAAli9fTvGldiIjIxk3bhzdunVDr9fj7OxM//79+fHHH6vzthT1BCUSlcDFzpLX/L0AWPRbNLfSsx58gkLxADIzM5kyZcpD7c6dO0fXrl1ZvHgxr732Gps3b+a7777Dy8uLkSNHMnz4cPLy7s/h2b59O1u2bOGFF17gv//9L8uXL8fDw4MhQ4YwZ86c6rwlRX1ASlkvDiBFr9fLmiYlPUv6TtkqPUN/ljO2nK3x6zcUYmJiZExMTNU1+NvnUl7aXXr9pd2aTQ2wZMkSCchBgwZJMzMzeebMmSL1Xbp0kQEBAVJKKXNycqSPj490cHCQFy9eLNHWrFmzJCCnT59eUHbz5k2Zl5dXwnbAgAHS2dn5of5V+WevqFXo9XoJpMhSnq0qkqgkehsLXu/fGoClf8Rw8/Y9E3ukKBPu3eH7URC9p2Rd9B6tzv3BO7tVNaGhoej1eiZPnlyqzfr16zlz5gxhYWF4e3uXqB83bhwdO3bk008/JTtby5O5uLgY3QOlZ8+eJCUlkZGRUXU3oah3KJGoAl7z98LRxoKM7Fy+3lVytzBFNZKbA7cul/9w8IRBM2Htq3Di+/vlJ77XygbN1GzK225uxYdDOzg4EBYWxoYNG9i3b59Rm/DwcACj+0yAtptiSEgIycnJHD58uNRrSSnZuXMnrVu3xtrausI+K+o/ame6KsDeyoK3AryZvuUc3x24zOuPtaKZXv3j1QhpcfBF58q18cNfy1ZWFsaeAMeK77swZswY5s6dy6RJk9i9e3eJ+itXtMmbrVq1KrUNLy+vAls/Pz+jNl988QWHDh1i8eLFFfZV0TCoUCQhhPiHEEIKIY4ZqQsSQuwXQmQIIW4IIeYLIUpsZiGEsBNCzBVCxBtsDwkhjH89qgO82scLFztLsnLymLfjoqndUdRRrKysmDp1Knv27GHTpk0VakMaRjeVts3uhg0bGD9+PKNGjeK1116rsK+KhkG5IwkhRCcgFCixj6cQYgCwGdgAvA80B2YCvkKI/lLKwsumrge6AxOBaGAUsF4I8YyUcnN5/TI11o3M+PtAb6b+dIa1h67yVoA3LZxsTO1W/aexu/btvTJcjbgfOTy3EFr0qpw/lWTkyJHMnj2bsLAwBg8eXKSuZUttpn90dDQdOnQwen7+TnItWrQoUbdp0yaGDRvGc889x8KFCyvtq6L+Uy6REELogEXAQuARoHiE8C/gFDAsXxCEEPHAduAFYI2h7EkgEHhOSrneULYTaA3MRhOaOsfwXi1ZsOcS8amZzP31Ap++0MXULtV/zMwr1b1D9B7YGgojf9Lefz8KXlgKrR6rCu8qhE6nY9q0aQwZMoQVK1YUqQsMDGTBggVs3LjRqEhIKfnxxx9xcnKie/eiifctW7bw3HPPMXjwYFasWFFkH2yFojTK2930DuABvFe8QgjhDvQElheOGKSU4UAcMLSQ+bNAKrCxkJ0ElgEdhBA+5fSrVmBlYcbbj7cFYN2RWC7dvGNijxQPJH8UU74otHpM+7m0UU81SEhICP7+/nz44YdkZd2ff/Pss8/SoUMHZsyYQVRUyUESn332GWfPnmX8+PE0atSooHzbtm08++yzBAYGsnbtWiwsLGrkPhR1nzKLhBCiNfARMEZKmWbExNfwespI3clC9fm2Z4p1PwGcKFRf/PopDzoAfVnvpTp54VEPWjrZkCdhzi8XTO2OojSKC0Q+tUgoZs6cSUxMDGfPni0oMzc3Z926dVhbW+Pn58e//vUvdu3axaZNmxg1ahQTJkxg6NChTJw4seCc33//nWeffRZ3d3cmTpzIkSNH2L9/f8Fx754atq0onTKJhNAyYN8C26SUG0oxcza8JhupSy5Un29bmh3FbOsUFmY6xj6hRRM/nbjGuevG9FRhcuKOlN6tlC8UcUdq2qsi+Pv7Gx3q6uPjw7Fjxxg5ciQLFy5k0KBBDB8+nKioKBYvXszatWuLdCX98ssvZGRkcOnSJQYMGECfPn2KHPHx8TV5W4o6hsgfCfFAIyHeAD4FfKSUcYayXYCDlLKr4f1LwArgUSnl4WLnrwCekFI2NbyPBM5LKZ8pZtcWiAT+V0r5TbluRIgUvV6vT0lJKc9p1UJuniT4891E3UznT53cmD/iUVO7VOfJT8Z6elYi/6CoEOqzr984ODiQmpqaKqUsMQoVyhBJCCFc0BLS04F0IYSDYUirOWBmeG8FJBlOMRYFOFE0ckh6gB0YjzLqDGY6wTtB7QDYdjqBk7GpJvZIoVAoKkZZups80Pr7pwO3Ch3+aLmDW8A/gNMG+xL5BLSRUIVzFaeBjobRUsXtwHheo07xpG8zOjS1B2B2+HkTe6NQKBQVoywicREYaOQ4DkQZfl4gpYwFDgEvF374CyGeANyBHwq1uR5t+GyR7ibgVbRuqDMVuptahE4neDe4PQC7zt/k8OU6HRwpFIoGykPnSUgp7wC7ipcbRhQhpSxcF4o2J2KVEGIB9yfTHQC+L2S3GdgJLBJCOKNNphsJ9AOGVOA+aiWBHZvQxUPP8dhUZm+PZOXrxpdIUCgUitpKlS7wJ6XcATwNeAGbgM8Mr4OllLmF7CTwZ2A1MA3YAnRGm1z3U1X6ZEqEEIwzRBN7o5LYezHRxB4pFApF+ajwAn9SygGllG8Ftpbh/DRgjOGotzzW1oWeXo4cjLnF7PBI+ng7l7qmjkKhUNQ21FLh1YwQ93MThy/fYlfkTRN7pFAoFGVHiUQN4NfamX5tXACYvf18iT2IFQqForaiRKKGGBeszZs4FZfGttMlFtBVKBSKWokSiRqie0tHnujQBIDPws+Tm6eiCYVCUftRIlGD5M/Cjky4w88nrpnYG4VCoXg4SiRqEF93PU8+0hTQVojNyS2+CK6iIbJ06VKEEFhbWxMbG1uivmvXrgwYMACA+Ph43nvvPfz8/HB2dkav19OzZ0+WLVtGXp76e1JUPUokaph3AtshBEQnprP+aJyp3WmwLD61mIj4iFLrI+IjWHyqZvd/zszMZMqUKQ+0OXz4MMuXLycwMJDvvvuO77//Hj8/P0aNGsW4ceNqyFNFQ0KJRA3T1s2eIV2aA/DFrxfIylHf/kyBr7Mv43ePNyoUEfERjN89Hl9nY8uQVR+DBg1i2bJlRfaPKI6/vz9RUVF8/PHHDB48mODgYL788ktGjRrFvHnzqA2rICvqFxWeTKeoOGMD2/HTiXhib2Ww9tBVXvFTSzBXlJy8HBLuln+0mLu9O5N6TWLcrnFM7j2ZLk20rWaP3zjOtAPTmNx7Mu727sTdKV+052bjhrmuYv9WoaGhREREMHnyZNavX2/UxtHR0Wh5z549Wbp0KfHx8Tg4GF3xWaGoEEokTEArF1ue7+7BmkNX+XLHBZ7v4YGVhdpvuCIk3E1g0LpBlWoj9LfQMpWVha1Dt+Ju516hcx0cHAgLC2PChAns27ePPn36lPncHTt2YGtri5eXV4WurVCUhupuMhFvP9EGCzNBQto9Vhy4Ymp3FLWEMWPG0KJFCyZNmlTmc9avX8+6det49913sba2rkbvFA0RFUmYCA9HG4b3asl/9l3m610XebFnC2wt1a+jvLjZuLF16EOXCnsgx28cL4gcZvafWdD1VFF/KoOVlRVTp05l9OjRbNq0iaeeeuqB9vv372fEiBEEBgbywQcfVOraCoUx1FPJhPx9YBvWHLxK4p0slu2L4W8D2pjapTqHuc68wt07oCWpZ0TMYFHwIgDG7x7PrIBZ9GrWq6pcLDcjR45k9uzZhIWFMXjw4FLtDh48yKBBg+jWrRsbN27E3Fz9OyuqHtXdZELcGlvxah8taT1/9yXSMrNN7FHDIn8UU74o9GrWi1kBs0od9VRT6HQ6pk2bxsmTJ1mxYoVRm8OHDxMcHEzHjh3ZvHkzNjY2NeyloqGgRMLEvBXgjU0jM1Izsln0W7Sp3WkwFBeIfGqLUISEhODv78+HH35IVlZWkbqjR48SFBSEt7c3W7duxd7e3kReKhoCSiRMjLOdJaP9WwGw+PdobqVnPeQMRVVwKulUqd1K+UJxKsm0W63PnDmTmJiYIvMmzp8/T1BQEDqdjo8++oizZ8+yf//+giMtLc2EHivqI6oTsxbwev/WLNsXw+3MHBb8donQQR1M7VK9Z7Tv6AfW53c/mRJ/f39CQkL48ccfC8r27dtHUlISgNGk9s6dOwuW8FAoqgJRX/Y2EEKk6PV6fV2dcTr31wt8Fh6JtYUZeyYOxNXe0tQu1SouX74MgKenmnhY06jPvn7j4OBAampqqpTS6CxM1d1US3jN3wtHGwsysnP5eleUqd1RKBQKQIlErcHeyoK3ArwB+O7AZeJTM0zskUKhUCiRqFW82scLFztLsnLymLfjoqndUSgUCiUStQnrRmaMGahFE2sOXuVq8l0Te6RQKBo6SiRqGcN7t6S53oqcPMkXv14wtTsKhaKBo0SilmFpbsbbT7QF4IcjsVy6ecfEHikUioaMEolayPM9PGjpZEOe1LY5VSgUiuLU1O6KSiRqIRZmOsYaoomfTlzj3HU1i1ahUBSlpnZXfKhICCH6CiG2CSHihBCZQoibQogdQogSy1MKIYKEEPuFEBlCiBtCiPlCiBITNIQQdkKIuUKIeIPtISFESKXvph7x527ueLvaIiV8Hh5pancUCkUto/A6Y3ti93DtzjWg9HXJKkpZluVwBM4DS4DrhvdvAJuFEMOllKsBhBADgM3ABuB9oDkwE/AVQvSXUhbezHk90B2YCEQDo4D1QohnpJSbK31X9QAzneCdoHaMWXmUbacTOBmbyiMeelO7pVAoTECezOPm3ZvE3onl6u2rxN6OJfZOLLG3Y8nJy+Hvv/6dprZN+cT/kypf7r5Cy3IIIczRHu4XpJSPG8oiAAugR74gCCGCgO3Ai1LKNYayJ4FNwHNSyvWGMgH8BjhLKTtW6Ebq+LIcxsjLkzw59zfOXb/NgPauLH3NtGsJmZL6vDTE0qVLee2117CysuLChQt4eHgUqe/atSsODg7s2rULACEEY8eOZc6cOSXamjNnDu+88w7R0dEFW5kuWrSIxYsXExkZSWpqKm5ubvTr148PPvgAHx+fh/pXnz/72kRGTgZxt+OMCkHcnTju5d4rUzuLgheVSyAetixHhRb4k1LmCCFSgWwAIYQ70BN4t3DEIKUMF0LEAUOBNYbiZ4FUYGMhOymEWAYsEEL4SCnPVMSv+oZOJ3g3uD2v/+cQu87f5FBMMo96OZnaLUU1kZmZyZQpU1i0aFGVtpuYmEhgYCATJkzA0dGR6OhoZs6cSe/evTl69Cht2qjNrmoCKSWJGYkFD/7CQnD19lUSMxLL1I6rtSse9h60sG+Bh50HHvYepGWlMSNiRrX4XWaREELo0HIYTYA3gXbAeEN1fnbE2NrKJwvV59ueKdb9BHCicH1Z/arvBHZsQhcPPcdjU5m9PZJVb/iZ2qV6QdLChVj5PoKtX2+j9en7D5B56iTOf/1rjfk0aNAgli1bxvjx4+nYsUIBtVFCQ0OLvA8ICMDPz4+OHTuycuVKPvzwwyq7VkMnMyeTa3euFY0GCkUEmbmZD23D0swSdzt3TQTsPQqEoIV9C5rbNcfavOg+5hHxEXx68NNq212xPJHEWrSIACAN+B8pZf7mws6G12Qj5yWj5R8oZGssE5tcqL4EQoiH9SPVyw57IbRo4tXFEey7lMTei4n0beNiarfqPFa+jxD3zju4f/55CaFI33+goK4mCQ0NJSIigsmTJ7N+/fpqvZaLi/Y3ZGFhUa3XqW0sPrUYX2ffUh+gEfERnEo6VepS8lJKkjKTikQAhYXgxt0bZfLD2cr5fjRQTAhcrF3QibINPDWWpM5PZtdk4jqfiWiJ6KbAS8BaIcRIKeWqQjalJTiKlz8oEVI/1i6vQvq3daGXlxMRMcnMDo+kj7czWhpHIXNyyL6eUO7zLDw8cHvvPWLHjqXpBx9g3bUrABnHjnH9n/+k6QcfYOHhQVZsXPnabeqGqOBe0w4ODoSFhTFhwgT27dtHnz59SrWVUpKTk1OiPC+veIB+n9zcXHJycoiJiSE0NBQ3NzdGjhxZIV/rKvnDRo09QPMfuNP7Tyc6NbqkEBiigYychy++aaGzwN3OvUS3UAv7FrjbuWNjUfntZsuyu2JVCEWF95MQQvwE+AMuQBCwFfiTlHJ7MbstgKeU0sfwfh9aGqJvMbvewH5gmJRybQX8qXeJ68IcuJTEsAX7AVjyWk8Gtm9iYo9qltKSp1mxcUQFBprCJaN4//ILjTzcy3VOfuL66NGjdOjQgXbt2tGqVSt2794NGE9cP4zCiet8XFxcCjYsateuHRs2bChTt1Z9S1znP1zHPzoec505sXdiOZJwhAPxB7C3tCclMwVZhu+qTlZOBQ//4tFAE5smZY4GKkplo6J8qiVxne8D8DTgCpw2lPmijWYqzCPA3kLvTwNDhRC6YnmJRwyvpt0zspbSu7Uz/dq48PvFRGZvP8+Adq4qmqiHWFlZMXXqVEaPHs2mTZuM7j4HMHz4cMaNG1eifNWqVXz22WdGz/n111/JyMjg0qVLzJkzh4EDB/Lrr7/SqVOnKr2H2kxUShSHbxzG2tya9/54r0T9rcxbBT+b68y1aKDQw7+wKNha2Nak6yWoqd0VKyQShiGrA4AUIMkw2ukQ8LIQYk6hIbBPAO7AD4VOXw/8BXiGQiOcgFeB82pkU+mMC27H7xcTORWXxrbTCQzybWpql0yORVM3vH/5pVJtZBw7xrXx2hiM5rNmFXQ9VdSfyjJy5Ehmz55NWFgYgweXmLMKQJMmTXj00UdLlP/++++lttulSxcA/Pz8CAkJoW3btkyePJmNGzeWek594OKti2y/vJ3tMduJSjW+oVdPt550bdK1ICLIjwbMdGY17G3t46EiIYRYAVwGDgOJQDNgJPA48LaUMr9jNBQtilglhFjA/cl0B4DvCzW5GdgJLBJCOKPNtxgJ9AOGVME91Vu6t3TkiQ5N+PXcDT4LP0+QjxtmuoYdTQhz83J37xQmff8BEj75hJZLl776HxkAACAASURBVAKUmsyuSXQ6HdOmTWPIkCGsWLGiWq5hZ2eHj48PkZH1bza/lJKLKfeF4VLqpSL1zWybEewZjLu9O9MOTAPgrS5vmXxP89pKWSKJfcDLaMNe9WhzHA4BIVLKn/KNpJQ7hBBPA1PRJsvdRpt9PVFKmVvITgoh/gxMMxwOaENenyvcnsI47wS149dzN4hMuMPPJ64xpGvFH5ANncKjmPJFwf3zz2uFUISEhODv78+HH36ItbX1w08oJ8nJyRw/fvyByfG6hJSSyFuRbL+8nfDL4USnRhepb27bnGCvYII9g/F18eXg9YOM3z2+2oaN1iceKhJSynnAvLI0ZhgSu7UMdmnAGMOhKAe+7nqefKQpm09eZ84vF3jqkWaYm6l1GsuLMYEAsPXrXWuEYubMmfTr1w/QupcqSteuXRkxYgTt27fH1taWyMhI5s6dy927d/nggw+qyt0aJ18YtsVsI/xyODFpMUXq3e3cCfYMJtgrmE7OnQpyeDUxbLQ+UZnEtcJEvBPYji2nrhOdmM4PR+P4n0dbmNqlOkfmqZOlikC+UGSeOmlSkfD39yckJIQff/yxUu34+fmxZMkSrly5QkZGBm5ubgQEBLBmzRp8fSu/SmhNIqXkXPK5gojhctrlIvUedh5axOAVjI+TT4nBHTU1bLQ+UeEhsLWN+j4EtjjvrDnG+qNxuDtYs3P8ABqZ1+9oor4Nw6xLmPqzl1JyNvks22M0Ybhy+0qR+hb2LQoiho5OHR846q+qho3WJ6pzCKzChIx9oi0/Hr9GXEoGaw9d5RU/9fBU1B+klJxJPsP2GC35HHsntki9Z2PPAmFo79i+zMPBa2rYaH1CiUQdxcvFlhd6eLD64FW+3HGB53t4YGWhhusp6i5SSk4nndaE4fJ24u4Une3u1diLIM8g/uT1J9o5tlPzhGoIJRJ1mDGPt2HdkVgS0u6x4sAV/tKvlaldUijKhZSSU4mnCnIMxoQhf1SSEgbToESiDuPhaMPwXi35z77LfL3rIi/2bIGtpfqVKmo3UkpOJJ4oyDHEp8cXqW+tb10gDG0c2ihhMDHqiVLH+fvANqw5eJXEO1ks2xfD3waovQEUtY88mceJmycKIobr6deL1Hvrve8Lg6P6G65NKJGo47g1tuLVPp58+1s083df4hU/TxpbNazlnxW1k3xhyJ/HkHC36Gq9bRzaFAiDt4O3ibxUPAwlEvUAawszLM11pGZks+i3aN4Jalekfm9UIidiU3krQP0jKqqXPJnHsRvHCiKG4vsrtHVsWzAqqbW+tYm8VJQHJRL1AD9vZxb8pq1Ps+j3aEb19cLRthGgCcSYlUeZ91I3U7qoqMNk5mSSlpWGJ8aHWe+P3094TDjmOnN+ufwLNzKKCkM7x3YFwtBKrwZX1DWUSNQD+nq78OWL3Xlj+SHu3Mth/p5LTBrcoYhA9PVWu9kpKoaZMONA/AGkgyyYQ5Cbl8uRG0dYfmY5u67uKrH/QgenDgR7BhPkGYSX3ssEXiuqCiUS9YSgTm4838OD7w/Hsvj3aNq52fHxprNKIBSVxsLMgt7NevPu7ncZ12Mcp5JO8euVX0nMSCxi19GpI8FemjB4NlaTO+sL9XsthwbGlJBONLYyJys3j3Frj/Pekx2UQNQBli5dihACa2trYmNjS9R37dqVAQMGABAfH897772Hn58fzs7O6PV6evbsybJly0rdunTTpk0MGjQIFxcXLC0tadWqFW+99VbBchtlwcXahcc8HuODvR+w5vyaAoHwbOzJ/3X/PzY9u4m1z6zlr4/8VQlEPUOJRD3CztKcCX9qX/D+nz+fJermHRN6pCgPmZmZTJky5YE2hw8fZvny5QQGBvLdd9/x/fff4+fnx6hRo4zuVDdhwgSefvppLCwsmD9/Ptu3b2fChAls3bqVLl268McffzzUr5y8HH698isbo4puTjS9/3R+fvZn/vLIX2jZuGX5blZRZ1AiUY/YG5XI579c4J9DOmFtoSMlI5vn/r2XizeUUBTnyLbLxJ6/VWp97PlbHNlW9m/aVcGgQYNYtmwZZ8+eLdXG39+fqKgoPv74YwYPHkxwcDBffvklo0aNYt68eRRe4HLVqlXMmjWLMWPG8NNPPzF06FACAgL429/+xsGDB3FycmLo0KGUtihmnszjevp1bmfdJvVeKjqhI9gzuKC+iXXD2me9oaJEop5QOEk9oo8Xa97sg42FNixWE4rbpnaxVtHEqzHbvj1lVChiz99i27enaOLVuEZ9Cg0NRa/XM3ny5FJtHB0dsbAoOQ+mZ8+e5ObmEh9/f/byJ598gpOTE//6179K2Lu6ujJ9+nQSEhJYtGhRifr07HSiUqJIykgCoHGjxoT1CuPg9YMsCl7EouBFjN89noj4iIrcqqIOoUSiHmBsFFNnDwfWvtUX20ZmpGVqQnEhof4JRV5uHmmJGeU+Gjtb0f9/2rJ1/kkiI64XlEdGXGfr/JP0/5+2NHa2Kne7ebnG8wJlwcHBgbCwMDZs2MC+ffvKde6OHTuwtbXFy8sL0HIXp0+fJigoqNSd7Z5++ml0Oh3h4eEFZbl5uVy7c42Y1BiycrMQCKzMreji2oV/H/t3wV4LhfdfUEJRv1Gjm+oBJ2JTjY5i8nXXs+bNPgxbsI+0zBxeXLCfla/70b6pvYk8rXru3LrH8vfL90AtTvjiM2UqKwsjPu5DY5eKbzc6ZswY5s6dy6RJk9i9e3eZzlm/fj3r1q0rstXplSvangutWpU+L8HW1hZXV9cC2ztZd7h25xrZedkAWJlb4W7nTmxaLBHXI9RGPQ0UFUnUA94K8C51FJOvu561b/bBwcaCpPQsXvp2P+eup9Wwh4qyYmVlxdSpU9mzZw+bNm16qP3+/fsZMWIEgYGBFdqKVEqJEIK423FcTrtMdl42QgjcbN1orW+NlbkVuTKX3s16GxWBfKE4lXSq3NdW1A1UJNEA6NRcz8q/+vHywv0GoTjAir/2pmOzmu1zrw7sHC0Z8XGfSrVx/VJqQeQQNNqHpq31lfKnsowcOZLZs2cTFhbG4MGDS7U7ePAggwYNolu3bmzcuBFz8/v/zi1baqONoqOjSz0/PT2dxMRE2vm2I+Welry2sbChuV1zLM3u34eVuRVW5laltqM26qnfKJFoIPg0b8zK1/14eeEBkg0RxYq/+uHTvG4Lhc5MV6nundjzt/ht7QWGvKMtW7Lt21P86XVfPNo7VpWL5Uan0zFt2jSGDBnCihUrjNocPnyY4OBgOnbsyObNm7GxsSlS36xZM3x8fNi+fTsZGRkl8hLZedksXbuUvLw8/AL80AkdbjZuOFo5qqW5FUVQ3U0NiI7NGrPqdT+cbRtx6242Ly3cz+lrqaZ2y2Tkj2LKFwWP9o786XXfUkc91SQhISH4+/vz4YcfkpWVVaTu6NGjBAUF4e3tzdatW7G3N55jmjx5Mrdu3SI0NLSgTEpJSmYKhy4eYsbUGTi5OjFi1Ai8HbxxsnZSAqEogYokGhjtm9qz6g0/Xvp2P4l37nc9+bpXvIulLlJcIPIpLBSmjihmzpxJv379AGjSRJuTcP78eYKCgtDpdHz00Ucl5lT4+PjQuLEWHb788sscPHiQL774gpiYGF559RWwhZOnTrJk3hJSb6Wy+ofVPNLiESUOilJRkUQDpJ2bPate98PFzpLUjGxeXniAk7ENK6K4EZNWqgjkC8WNGNMm+P39/QkJCSlStm/fPpKSkkhKSuKpp56iT58+RY4jR44UsZ8zZw4bNmzgTsYd3njjDV4e8jKL5i4i4IkADh89zFOBTymBUDwQIaV8uFUdQAiRotfr9aXNHlWU5OKN27y44ACJd+7R2Mqc7/7am84eDqZ2yyj56wx5eqp1gcrDvZx7XEu/xt3suwCY68xpZtuMxpZlz0Wpz75+4+DgQGpqaqqU0ug/v4okGjBtmtiz+g0/XO0tScvM4eWFBzh+VYlsfUBKSWJGIlGpUQUC4WDpgLeDd7kEQqFQItHAadPEjtVv+NHE3pLbmTm8svAAR6+YNmmrqByZOZlcSr1EQnoCUkosdBZ4NvbE3d4dc51KQyrKx0NFQgjxhBBiqRDivBDirhAiVgjxgxDiESO2QUKI/UKIDCHEDSHEfCFEiRBGCGEnhJgrhIg32B4SQoQUt1PUDN6umlC4Nbbk9r0cXl0UwRElFHWOPJnHjbs3uJRyicycTACcrJzwdvDGrpGdib1T1FXKEkm8BbQEPgcGA+MM7w8KIfzyjYQQA4DNwFXgGWA8EAJsEkIUv8564GXgfeAp4AywXgjxZGVuRlFxWrvasfqNPjRtbFUgFIcvK6GoK9zNvsullEvcvHsTiaSRWSO89F40s2uGmc7M1O4p6jAPTVwLIZpIKW8UK3MAooEdUsqhhrIIwALoIaXMM5QFAduBF6WUawxlTwKbgOeklOsNZQL4DXCWUnas0I2oxHWVcDkpneEL9nMtNRPbRmYsG92LR72cTO2WSp6WQn70kL9aK2gbBLnauKIr8d2sYqjPvn5T6cR1cYEwlKUAFwAPACGEO9ATWJ4vEAa7cCAOGFro9GeBVGBjITsJLAM6CCF8ynBfimrC09mW1W/0wd3BmvSsXF5dHEFEdLKp3UKn05Gbm2tqN2oV6VnpXEy5WCAQluaWtHZojZutW5UJBEBubi46nUpfNlQq9JsXQrgCvkD+ql6+hldjq3ydLFSfb3umsJgYOFGsreLXTHnQATSs2WDVSEtnG1a/4Ye7gzV3s3IZtSSCA5eSHn5iNWJlZcW9e/dITja9YJmaguW802LIztUW5Gti04TW+tZYm1d8iRJjJCcnc+/ePaysSl+7SVG/KfdQB0PX0AI0gZllKHY2vBr7D04Guhd67wxElmJXuC2FCWnhpAnFiwv2E5eSwaglB1nyWk/8Wpvm1+Pi4sK9e/dISEggJSUFM7OG2c+enZtNek46+d3EZsIMWwtb7t65y13uVum1cnNzuXfvHvb29ri4qL3SGyoVGQ/3KfBn4DUpZfF9FktLcBQvf1AixGhdaf1l+ahooupp4WTDmjc1oYi9lcFrSw6yaNSjpS5LXp0IIXB3dycxMZHMzEzy8iq+uU9d5F7uPY7fPM6VNG3vBzNhRieXTrR1aFttM6YtLCwKBELNym64lEskhBCfAO8CY6WUSwtV5fdFGPua6UTRCCPpAXZgPBpRmAgPRxvWvNmHFxfs42pyBqOXHmTxyJ70bWMaoXB1da3x65oSKSXbYrYx/eB0kjO1f43eTXszpe8UWti3MLF3ioZAmXMSQoiPgMnARCnl3GLVpw2vxvIJj1A0V3Ea6GhkWGz+vAu1e0ktw93BmjVv9KGlkw2Z2Xm8tvQgv19INLVb9Z4bd28wdudYJuyZQHJmMnYWdkzpM4Vvg79VAqGoMcokEkKIKcAHwAdSyk+L10spY4FDwMuFH/5CiCcAd+CHQubrAQe0uRSFeRU4L6Ws2L6RimqluYM1a970w9PZhns5efxl2UF+u3DT1G7VS6SU/HDhB/684c/svLoTgACPANYPWc/z7Z5XXT+KGqUs8yTeRUtQ/wx8Uqz6npTyqMHucbQ5EevQEtvNgZnAFcBfSplrsBPAr0BnYCLafIuRaCIxREr5U4VuRM2TqBGup2Yy/Nv9RCem08hcx7evPkpAu4bVBVSdxN6OZeq+qeyP3w+Ao6Ujk3pNYnCrwUocFNXCw+ZJlEUkdgEBpVRfllJ6FbIdBEwFugC3gQ1o3VNFpu4KIRoD04Dn0aKKM8BHUsoNZbin0vxUIlFDXE/N5KVv93PJIBQLRvRgQPsmpnarTpObl8vq86v54sgXZORkADDYazCTek/Cycr0kxkV9ZdKi0RdQYlEzZKQpkUUl26m08hMx/wRPRjYQQlFRbiUcokpe6dw7OYxAJpYN+F9v/cZ2HKgiT1TNATUUuGKasGtsRWrX/fD29WWrNw83lx+mB3nEkztVp0iOy+bb098y/M/PV8gEEPbDmX9n9crgVDUGlQkoagUN25n8tK3B7h44w4WZoKvX+5BoI+bqd2q9ZxNOsuHez/kXPI5ANzt3PlH33/g18zvIWcqFFWLiiQU1UoTeytWve5H2yZ2ZOdK/nfFYcLPqIgCYPGpxUTERxQpu5d7jy+OfMHwTcMLBOKVjq/wQ8gPSiAUtRIlEopK42pvyao3/GjvZk92ruRvKw6z7fR1U7tlcnydfRm/e3yBUBy7cYwXfnqBhScXkitz0QkdYb3CCO0Vio2FjYm9VSiMo7qbFFVG0p17vLzwAOeu38ZcJ5j3UncG+TY1tVsmJSI+gnd3v0sPtx7suLIDiUSHjkZmjZgzcA7+7v6mdlHRwFHdTYoaw9nOkpWv+9GhqT05eZIxK4+w5WS8qd0yKbaNbLE0s+TXK78ikbS0b4ltI1u+euIrJRCKOoESCUWV4mTbiFWv++HTrLEmFKuOsulEwxOK3LxcFp1cxCubXiHhbkLB/g5Xbl9hzoA59GrWy8QeKhRlQ4mEospxtG3Eir/2plPzxuTmSf7f6qP8dPyaqd2qMa6nX+f18NeZc2QOOTKHVvpWvN/7fVO7pVBUCCUSimohXyh83TWhGLv6KBuPxZnarWpne8x2hv44lIPXDwIwrP0wJjw6gS+Pfsmi4EUsCl5UJJmtUNR2VOJaUa2k3s3mlUUHOBmXik7A58O6MqSru6ndqnLSs9OZETGDDRe1lWUcLR2Z2ncqtha2jN89nlkBswq6mCLiI0qUKRSmQiWuFSZFb2PBd3/tTRcPPXkS3llzjPVHY03tVpVy8uZJXvjphQKB6Nu8L+tC1hkVCIBezXoxK2CWiigUdQIVSShqhNSMbF5dHMHxqykIAbOe78LQHh6mdqtS5OblsujUIv597N/kylwsdBaM6zGOlzq+hE7oWHxqMb7OvqVGCxHxEZxKOsVo39E17LlCcR+1wJ+i1pCWmc2riyI4ZhCKT5/vwvN1VCiu3blG2G9hHLlxBIA2Dm2Y0X8G7Z3am9gzhaJ8qO4mRa2hsZUFy//Si+4tHZASJvz3OGsPXTW1W+VmS/QWnv/x+QKBeKnDS6x6apUSCEW9pFx7XCsUlcXeyoJlo3sxaslBDl++Rei6E0gpGdazpaldeyh3su4w7cA0frqk7YvlZOXEP/3/yWMej5nYM4Wi+lCRhKLGyReKRz0dkRJC151kVcQVrfL3ORC9p/STo/doNjXMsRvHeP6n5wsEor97f9aFrFMCoaj3KJFQmAQ7S3OWju5FLy9t17WwH06y8sAVcO8O348yLhTRe7Q69+415mdOXg7/PvZvRm4dSdydOCzNLJncezJfPfEVLtYuNeaHQmEqlEgoTIadpTlLXutJr1aaUExef5LvEjzhhaWaGFzafd84XyBeWAqtaubb+9XbVxm1dRRfH/+aPJlHO8d2rH5qNcM7DFf7TSsaDGp0k8Lk3M3K4Y0lf3Ar5iS+uhhGt06lfcYxSIwEK0do1gXiDsEzc+GRodXuj5SSny/9zCcHPiE9Ox2AET4jGNt9LJZmltV+fYWiJlFDYBW1j6x0SDgN8cch/hjEn0DeOIvIy37oqRn2nlyyf5ROfZ+GVgFg61ylrqVlpfHx/o/ZEr0FABdrFz7x/4S+7n2r9DoKRW1BiYTCtGSmwvWTBkEwHImRIPOMmt8SDhzN8eKU9KJfu6Z0v/SNVmFhC4Zv9UVwewRaB2iC4dkXLO0q7OrhhMOE/RZGfLq2au2AFgP4qO9HOFo5VrhNhaK2o0RCUXOkJxYVg/jjcCu6dHt9S2jWWetOMhyZVq688M1e7OL3Mc9iLts6ziCgvStNtr7Jp1lDebmLA56pB+HqAcjJLNqezhzcH70vGh49wbzRQ93Ozsvm62Nfs+jUIvJkHlZmVkzoOYEX2r2gcg+Keo8SCUXVIyXcvl5SENIesCaTk3cRMaBZF7BxMmqadWEXGStH8Na9t9mX1wmAvrrTzLeeR1zg17Tt/SRmufc0oYjeA9G7Ie4IyNyiDVnYQEs/TTBaB0DTzqAzK2JyJe0Kk36bxMnEkwB0dOrIjMdm0FrfuuKfj0JRh1AioagcUkLK5ZKCkH7TuL3QgWuH+0LQtDM0fQSsGpfteoZRTFnPLuadCHs2nby/V3Yf3WnmWcwlzOxd7Ds8TpBPE/q3dcXW0lzr1rq8VxsRFb0bbpwp2baVA7TqD60CkK0C2JBymukR08nIyUAgGOU7ire7vo2FmUUFPiiFom6iREJRdvJyISkKrp8wJJQNgpCZatxeZwFuPoWig67QxAca2VTch9/naPMgWj3G3qhEXvr2AADPdmvOuet3cEjYR2dxifm5zwDQyFyHv7czgT5uBHZ0w62xldbOnRv3o4xLuzWhM5Cq0zHVxYlwW83PJhb2TOv1Hr3bPFVxvxWKOooSiYZAoQerUaL3aN0x/f7vflluNtw8XzQ6uH7SeHIYwNxaiwgKdxe5dihTn39F2BuVyJiVR5n3UjeAgp9bOtnw69kbhJ9JYP+lJHLyiv79dvbQE9jRjSAfNzo0tb+fU7gVA5d2E3HxJyannyXBTJsiFJh+lymJyTjk5YFzm/tdU179S+0OUyjqE1UiEkIID2AC0APoCtgCA6WUu4zYvgSEAu2BROA74B9Sysxidm7Av4CnAGvgCBAqpdxb5rsr2l7DFYkHTTTLrxs4GYTZfUFIOA2594y3Z9lY6yYqLAgubUv051cXhQWir7dLqWVpmdnsPn+TX84msPPcDdIyc4q04+5gTZAhwujW0p4Fp/7NklNLkEiszSwJc+3LnxMTEJf/gKzbxbwQWlI9XzRa9oFGtjVx+wpFjVJVIjEAWIv2IL8HhGBEJIQQrwDLga+B74GOwExgk5TyxUJ2VsAhwA4IA5KA/wMGAn2llEfLdZc0cJGA+2Lw3EKtuyf+OERug0u7AFnqkFNsnEsmlB28QGeayfjGxKAsddm5eRyMTib8bALhZxKIvZVRUCca3cTOYzVYatundnTqxKcBM/Fs7KkZ5ObAtaMQvUvrmroaUVJAdRbaaKmCkVOPgspdKOoBVSUSOim1p4wQ4s/AeoqJhBDCDIgFIqSUQwqVvw4sAPyklAcMZX8DvgJ6SCmPGMosgbPAeSnl4PLeaIMXCYCj38GPb5cuCPbNDUJQKEpo7A61aJjnN7uj6OyhLyEC+eyNSuREbCpvBXiX2oaUksiEO4Sfuc66C+u42WgtQpeNlIKspAHkJQXRq5VrQbdUC6diOZTsDLiy/35O49rRkp+pha02LyNfNNx8SwprRboBFYoapspzEg8QCX/gd2ColPKHQuU2QAowR0o50VAWDrhJKTsXa/sTtK4qRyll8fj/YX41bJG4EA7/HQ330u6XtewLbZ7QEsrNOoNdE9P5V8PcyrzFP/b+gx1XdwDQ2NyV5tmjOXHRmXs5RR/4HZraE9jRjUAfNzq769HpiolmRgpc/uP+yKmb50pe0NqpYOQUrQeAU2uI+e3h3YA1uBaVQmGMmhSJN4FvgPZSyshi55wFoqWUTxrexwM7pZQvFbMbBqwGekspy7X5b4MVCSnhwDewbbL2bdfaGTKStLqRPzXIB9Dea3t5//f3uZmhDdMd5DWI9/3eR2+p525WDr9dSOSXMwnsOHeDpPSsIue62lsS2LEJgR3d8G/jgpWFkTzM7euFRk7tgdQrJW0ae2hRhp0bHF4C//Of+78LJRCKWkRNisRk4BPAVUqZWOycPwBzKWVvw/ssYIGUckwxuyBgO/CklHJLsbqHPf31er2eBiUSudmweTwcXqq9d24Dd5O0BxI0uAdRVm4WXxz5gv+c0e7fxtyG9/ze45nWzxidOZ2bJzl29RbhZ27wy9kELt64U6Te2sKM/m1dCPRx4/EOTXCxM7K4n5TarPL8KCN6j/Y7KI4wg64vQfsn4ccxDer3oqjdPEwkqmNnutJUp3j5g9SpfozLrU7uJsPaV7UuDQCvfpBwpug31vwltxvAAykqJYrQPaGcv3UegM6unZnRbwYtGrco9RwznaCHpxM9PJ2YNLgD0Ynp/HImgfCzCRyKSSYjO5ftZxLYfiYBIaB7S8eCPIa3q60mPEJoXUtOreHR1yAvD26cNojGHq2bKuuONhv86HLt6PQcuLSrqY9GoagUqrupLnIzElYNg+RL2vvOw+Hi9gbZ9y2lZM35Ncw6NIt7uffQCR1vdH6DNzu/ibmu4t+BbqVnsfO8Nh9jd+RN7mYVXfKjlYttQbdUD09HzM1KGQ2Wm60lp/d9BWc33i83a6RFFn3/HziXnoRXKKqbmuxu6gf8RtkS178ATYwkrj8GJgFOUso0ykGDEYmoHbB2FNxLBXMrGPIVpMY2yFE0SRlJTNk7hd2x2uZE7nbuTO8/nW5NulXpdTKzc9l/KYlfzibwy5kbXE8rurCgg40Fj7dvQpCPG/3buWJnqYlTwUgt3RlNqJ/5Es5vhmMryA+W89Ch8wnRfjfNq9ZvhaIs1KRImANXgf1SymcLlf8FWAj0kVLuN5T9HZgHdJNSHjOUNUIbAntBSjmoXE7RQEQi4lvYEqp1Xdi5wYurwKOHqb0yCb/H/c77v79PUqbW//9066eZ3Hsy9o3sq/W6UkpOxaURfjaBX84kcCa+6HeZRmY6+hiWCXG0seDHDWv4ymIuFi8uuy/ikeHkrn2VWzkWuFBoyZPWA6HfO5pdLRqWrKjfVJlICCGeN/zYE5gI/AM4DaTnJ5mFECOBpWhzIP7L/cl026SULxRqywptYp4V2mS6ZGAs8ATQT0p5uFx3ST0Xidwc2BoKBxdq75t2huGrQe9uWr9MwL3ce3x++HNWnF0BgJ2FHe/7vc9TrU2z7lJcSga/Gibw7b+URHbu/f+nPrrTfN3oS8bmjuVPT7/A8F4tEUKwNyqRpSuW85XFF1j0Gg3ntxRdkLB5d00sOjxVY7PcFQ2XqhSJ0gwvSym9Ctm9Q024ogAAIABJREFUgjbXoR3ashwrgClSyozCJwkhmgKfoi3LkS8ak6SUv5fJoZL+1U+RyLhl2O95l/a+4zPw7PwGuURE5K1IQveEcjHlIgDdmnRjev/puNvVDrFMy8xmT+TNguG1w7PXc0K2LljuvIWjNcGdmvLDkVi+erm71g0VdwT8x8KF7fD753Bl3/0GndtodZ2HgbnaNlVRPagF/uoySVGwchgkXdDe938XBr5vsiUzaoLFpxbj6+xLr2a9CsqklKw8t5LPDn1GVl4WAsHfu/6dvzzyl0olp6uT7Nw8DsXcIvxMAuFnr3M1+f53JGsLHS/19mSEnydeLsXE/sp+baZ2ZKER4PbNoM/foccosKze7jRFw0OJRF0leg+sGQGZKdpImJB50GWYqb2qdiLiIxi/ezyzAmbRq1kvEjMSef+P9/kj7g8AdOiY1HsSwzsMN7GnZUdKyXf7L/PBxtNFyoWAAe1cebWvFwFtXYvO9E44A398ASe/v7+ZkpUeer0Bvd8CW+PLligU5UWJRF3k0BJtklxeDti6wosroUWvh59XT8gXipc7vszKcytJzkwGoJGuEZ8P/JzHPOrWUN7CCxOmZWTzzprj2DTSkZSeXWDj6WzDCD9PXni0BXrrQgsHplyBvfPgyH8gxxCNmFtDt1eg79vg6FnDd6OobyiRqEvk5sD29+HA19p7N18YvgocWprWrxokKzeLrTFbmX98Pldua8tdWJtboxM65g6cW6Qbqi5Q2rLnf19xhFF9W3EgOom9UfdnaFtbmPHnbu6M7OtJh6aFdvNLT4QD8yFigRZdgjaL23eoNnzWrVNN3paiHqFEoq6QmQr//QtcDNfet38SnvsWLO1M61cNkZyZzNrza1lzfs3/b++8w+Mq7r3/mbNdqy65yHKNbYyMDTHFxg7cEAimk9BCsQ2EQAKEwCVPbiA3ed+QRkhyLw4JaQQngAsQioFAEkJ9jXGh2OCKK5IsWW7q0mrrmfePObtarXZVbEm7kubzPOc558yZnZ0zWv2+U3/DkbYjnZ4vmb9kSAhEsmcjsl08sbaC5zZUdVi0N2dSITfMm8i500fhiC7WC7TAhsdV66J5f3uCU89TM6ImzB2IV9MMIbRIDAbq9sKKa+CIcinB5+6Cc344LKY/7qzfyfLty3l5z8sETeVszyZszJ8wn1mjZnH/+vuBwSkSvXV73uQP8fyHVTyxtoK9R9p3CByd62bBnPFcO2d8u/+ocBA2/00NckcnNgCMm6PEYup5Q3qCw9Gy4dUKRk7MZey0gqTPq3bUc6i8iZPPy/xuvL56Fy0SmU75u/D0QmirUxvbXPIQzFqQ7lz1K6Y0WV29mqXblrKuZl0sPNeZy5XHXcm1x19LZVNlbAAb6DCYPdQxTcnq3Ud4Ym05b3xyiOi/qNNmcNGJJVw/dwKfHZevfEeZJux4RU2frY5bXjSiTHVDzbhCb44UR9WOel798xbOu2VGJ+Pa1bNMpK/eRYtEJrNhKbx8N5ghtUPc1cuHdHeBL+TjpT0vsXz7csqbymPhE3MnsrBsIZdMvoQsR1anGU7QedbTcGFfnY9l6yp46v19NLa1D3SfODaP6+dO5OITS5Q7cymhfLUSiz1vtCeQN04NcM9apHYsPAqGUu0bkhvQvhQIaUpMKZERiWlKdW9KzEj7tbTu1TUJ9yamSezz8fE7fN6U1O5vZfvq/UybO5ri0mxOOLO01++iRSITMSPw+g9hzW/V/YgyuO4pKJiY1mz1FwdaD7DikxU8u/NZmuP2kp5bMpeF0xdyRukZGEJ1jXQlBsNVKADaghFe+riax9ZUsD3OFUih18k1p41jwekTKM33qMCaj1U31LYX2nfUyyqC2d+A2bdAVmGvvjvTat9mxCQcMgkHTcLBiHWtzpGgSSgYIRIyCYciVpy465CK13S4jf27GxgxPgdpwpF9zRSOycLpdsQZ6yQGPSKRUnaOE2fQ0+XD2pVl5/xvzOz130OLRKYRaIbnboad/1L3U+fDFUvAndv15wYhHx/+mGXblvFaxWtErLn+TsPJJZMvYUHZAqYWTO30mWSL6eJ5r+Y9ttRu4aYZN/Vr3jMVKSUfVtTz+NoK/rm5hrCp/n8NAedOH8UNcycyd3KR6oqq3aMqIh+taN+z2+FVi/LmfrNXbl26rH3ffAKjJ+d1MMiRkGWsg2YHIx4OdjTW4aBJJMHQd4wXLwAqzDSHhs3qDsMQCEMgbALDENa9FR4NsxkIAZGwSdMR5XjyS3fP6pVga5HIJOor4Mlr2v30nP5NmP+TITVAHTJDvFHxBku3LWXTkU2x8GJPMddMu4arpl1Fobt3NVlNcg41+Vm+vpIV71VyuDkQC586Mpvr503k8lmleF12aD6oplW/v6R9e1vDodx9fO5OGDGtQ7pmxMTXFKSlPkBrQ4CW+gAtDQEOVzZTs6sBT56TSMjE3xLCsAvMSPpqz4kIQ2B3GtgdBnaHTV07bdgdBjaHdW0997eGKd+kZtJNPXUUeSM9CCNqfAVCWOfYPRiWcRZxz6PG3Ohk0OM+nzJN0SFNQ0TTJulGWamo2lHPi4s3AlokUpLxIlG5Dp5aAL4jYNjhov9VNbohQmOgked2PceK7Ss46DsYCy8rLGPR9EWcP/F8HHoAtV8Ihk3+tfUAT6wp54OK+lh4jsvOFaeM5fq5E/jMiGzwNxJe/1da332alhZBa6SQlkgxLXmn0JpzEi3+LFrr/fiagvSlWbA5oka73UjbHDYc1lkZcsuod4hj4HDaLOMeZ/QTzrF4TgNbqn09EohvGQGDasA6kWN9Fy0SmcDHT8FL34JIEDwF8JWlMOnMdOeqT/i08VOWb1/OS3teos1aESwQnD3+bBaWLeSUUaf0qkakOXqCbWE27jjCP9ZX8fGuWjxhyDEF2aZgtNNBrimI+CPdJxSH02Mnu8BFdr4Lb76LSMRk53pVCTj1oomMmpjbwbDHG3Cb08BuNxBGZv39+3vgeiDpi3fRIpFOTBPe/AmsflDdFx+nXHwP8p3IpJSsq1nH0m1Leaf6nVi41+HlsimXcV3ZdYzLSb1t6HChr2YFSSnxt4Tau38aArTU+2NdQdGwUC8FwJPjINvdhjewm+zQp3httWTbavEW55F9+mV4T74Qp7fd++xQqH1n2iD8saCnwPaSjBOJQAus/AZ88rK6n3w2XPlX8CT9OwwK/GE/r+x9hWXbl8XcdYPaEW5B2QIum3IZ2c7hsUK8J/Tkn3j+16ZTMNrbqf+/NV4IGgKY4Z7/nwpD4M1zkl3gIivPRb2M8FFtC5tqW2g2JM1CEnYJLp1VyvVzJzK9JEfteLh6cfue6QAFk9T02bZ6qsxTefVlW/Ia68URxto/GhQ7Hw6l6bx6MV0vySiRaKxSA9QHNqv72V+H834Otsx0a90dh32HeWrHUzyz4xnqA+193qeMOoVFZYs4a9xZ2IbQ4HtfUrH1CK8t2caMz5fidNtpaQhwZF8zB/Y24fLYCPjCver/tzmMWNdPdoE6vPkusvPdeK1uIU+us6NHWYvdh1pYtq6CZz+soiUQjoWfNrGA6+dO5PwZo3HUbFBi8ckrREejq+TpvHr4Ns673M3YL3bcNLLq7bd59W+NnPeVPMaeddbRFJEmzWiRGGiqPoAnr4XWQ8oB2wW/UHPTByHbarexbNsy/ln+T8KmMip2w84FEy9g4fSFTC+anuYcZg5mxKTxcBu11a3U7W+hbn8rtftbaTzchuzhlE2nx95u/OOEoD3MjctrP+YxnpZAmJUbqnh8bQW7D7XEwkfmuFgwZwLXzhnHSH8lrHkIPn6aDU0XMdKxm7GurXDC5XD+/ZAzWrmzf+ZGqmb/hUOBSYOi9q3pjBaJgWTzs/DC7WpOujsPrnocJn8hffk5CiJmhLf3vc3S7Uv58GC7m4cCVwFXTbuKa6Zdw4isEWnMYXqRpqS5zk/tfiUGShRaqT/Y2mWXkMNlEAqohW0TTyxm1MRcJQBxguB0D2xLU0rJ2j21PL62nNe2HSSqZQ6b4IIZJdwwbwIn5/vY9OzPmVHzPLawz/qkAYWToKka5n6TD4su4v2GPG49a8qA5l/TN2iRGAhME97+Oaz6pbovnAzXPQ3FnReLZSotwRZW7l7J8u3LqW6pjoVPyZ/CoumLuHDShbjt7jTmcGCRUtLaEKSupr1VUFfdQt0BH+FA6gFiu8tGYYmXojFeCsd4KRqTTWGpl/qaVl789UdA7+exDwRV9T6Wr6/kqfcqqfe1u/+YUZrLvMnFvPr+NpaeuInx25dA3Kr5KEF3Ec4Jc2DsqTB2NoyZlZEejHvrdHE4oEWivwn64IVbYduL6n7Sf6gWRC9dH6SLfc37WLF9BSt3r6Q11O559MzSM1k0fRGnl5w+5KewtrUEqdvf2i4GVndRwBdO+Rmb3SB/dFZHMRjjJafQ3WnK52CaFeQPRXh5Uw2Prylnc3VjLNzrshExJX+d18Dc9bcC0FJwAtTvJZu2zgkJQ+1xMfY0JRpjT1Oz+tL8W+qp+/ZUIpJJ1D76KO4ZM/GePifp89Z16/Fv2UzRzTd3mY4Wif6kab8af6hRNURO+Spc+KuM97oppeTDgx+ybPsy3tr3Fqbl38dj93Dp5EtZULaASXmT0pzLvifYFqauJioG7S2EtqZgys8IQ5A/0kPhGC+FY7JjopA3woPRg4Vbg3VOvpSSjfsaeGJNOa9sriEUkcw1tvKw4zf8T+738Lps3H74p/x55A9wF5QwqW0r43xbGde6hWJ/edI0A448juSfSF3hZ2ksPInmopPAlYPdZmA3BHabwG4Y1rn92mET2IyEONa1w2ZgM1T8nlRmah99lD2F47l1q+i0EdQdKzbyxxMkk+squzWsmUDruvVU3303pYsXdxKKrp4lokWiv6jeoASi5YCqNZ3/gJrFlMG17lAkxL/K/8XSbUvZXrc9Fj4qaxTXlV3HFVOvIM+Vl8Yc9s20vnAwQv0BX0wIoqLQUhdI+RmA3GI3hVaLoMgShYJRWdgcR7cvw1CZk3+4OcA7/36eL2z+LrcHv8VaU+2CFxWNO0J3xsIAcmnls8ZuZondnGzs4rPGbvKEr1O6phTslGPZYE5ho5zKBnMqe2UJkqMr76hYRIVDiYsSleh12YFd3PDaIyw552Zed43jhNJcirxO1u6p49vFTcxb/iDOnzzA2C+cicfZvzP2pGkiQyFkMKjO8dfJwjo8DyFDQQI7d9GwciW5F16Aa8pUir56Y68EArRI9A9bV8LK29Sew65ctf5h6hf7/3tT0J1TvNcrXueZHc+ws2Fnh13fTiw+kUXTF3HOhHNwGJnR+umNYY1ETBoO+mJCEBWDpsNtXU4r9ea7KBrjpSBeDEZn9fnA8VCZk1/7y+/jrnkK1x1/4bd7S3joDbXJ0ZlTi5nDVq7e8GNWB+fzyhm3EIpIIqYkFDEJm5KwKYmEw5SEq5ga2k5Z+BPKIjv5jKzESOLwqVFm8ZE5hY1yChvMqXxkTqaJYxzbkBKHGcYVCeGMhDjpyG5u2/QCj5WdT0VuCQ4ZYWp9JVfvfJMXJp/Jfu8I7GaYXDsUOqHQaZDngDy7JMcG2TaJV0g8hokRDqUw5tHwYPt1grEn0rvFj91h5OUy9qHf9EogQItE3yIlrPoVvPUzdV8wEa59GkYe33/f2QNSudDeVb+LxR8u7rAq2iZsnDvhXBZOX8hJI05KR3a7JVEMTFOy870DrHpqJ5NPHkEkaFK7v5WGgz7lXC4Fbq+DolJvrHVQOMZLYYkXtzczBHGw0Prof1H9yNs03PuzWDcNEOueyX/g+5R+/Sy8N/+qR+lJ00Q2HsYsfx9Z/j5m5UfI6s2YviZkRGBGBDIikGF1HXaNJJw1hpBrFEF7EWGRhQwEMf1+ZCBgHX4IBGKHCEaPICIURAxyOyecTnU4HOpIuDb9foK71QLX8Y891mOBAC0SfUeoDV78Jmx5Tt1P+JzyweQt6p/v6yVRofjl539JMBJk2bZlrK1ZG3ue48xRu75Nu5aS7JI05jQ5QX+YlroAzXV+muv81OxqYM/GQ2QXummu9XcpBg63LdYiKIwbSPbkONI66N5XA4sDiTRNTF8bpq8Vs7UVs9WH2drKjjdWw4on4OIvM2HaRGTAT9VH2+GdtxCnzmbEmJHIgB/TH0D6/ZiB6NmP9AdiZ+n3I0Oh7jMywAQNG3aPB4fLiWl3ELHZCAs7QcNGQNjwY9AmDVpNgU8ahISNsGEnZNgJGzZChp2QYYtdx4e5PS6ysz3k5GSRl6uOgrwsCvOzKcr3UlyQjcOVQgScTrDZkv6O42dqta5bT+WNNwLtItHTmVpaJPqC5gPw1HXt20POWggXLQa7s++/qwdIKakP1FPZVElFUwUVTRVUNleyvXY7lc2VHeKOyhrFzTNv5tLJl5LlOLqdyY45v6bE1xSMCUBznb+DILTU+bucSRTF5jAoLPF2aBUUlWaTXeDKyBlYfTWw2BUxo95qGXWfr/062X13cdra6FMXsEeBcLkQbjeG04lwCAwRRuDHiLQgRBDDJhHWoa7ByM5HFI3DGDERMXIKxohJCE8Wwu3CcLsRLuu89Wl2uSZz+7ps/memnZz7vgtA8/0PcetWwRNn+5kh93TpYiQYNjnY5Kem0U9NY5s6N7Sxv9HPASvsSEvqyRCd3lfAiGwXJfkeSnLdlOS7GZPnoSTfTUmem5I8DyNzXNgTJkrED7bn/+o+ShcvBqD67rtp+K/7Og3Op0KLxLFS87EaoG6qBgTM/6nasGUAjFKDv4GK5goqmyqpbFaCUNlUSWVTJc2hznPVE7lr1l3cNPOm2K5v/UUoEKGl3hKAWj8t9YGY8W+uU/ddtQTiMWyC7AIXOYVuDJtg33blBuScG8s4bvbopO4mMplEMZCmScuqVey/515G3XsvrilTjsqYR3yqli99nQeE+wPpcmHPzgabjcihQwC4jj8e+6iR1IcF9WGDaeOLER43hssdZ5zdGG4XokNYwjO3G8MSBuHqQvClhPpy5dWg6j2oel+5vjGTVDAcXig92ZqCax3ZI9jy7t8pfe12qsd9G+ejz3QwrMGbr6J034NUn/t7ZnzukmMqL38o0kFI9jeo84FGf+w6fj1KdxgCRubECUiem9F5bhybNnD8Iw/Qcs+POOuaC7EZgvee/zfyR/+N+OH9zL58frdpa5E4Fra/DM/fAiEfOLPVDnLTzu/+c72gKdgUaxFUNlXGRKGiqYKmYFO3ny/2FDP7wAVkl9ooOS6PYDjIHzb9AYAl85cwu2T2MQ2QSlPiaw52qvlHjX9zrR9/a89/7K4sOzlFbrIL3OQUuskuVIIQPdqeWYpn5kzqC45LuragoH7ngHTRSNNU3SU+H2Zbm2WgfZhtPkyfDxkN80XP6pmMD7PihOvqiBw+DE6n6jMfAITHg+H1Yniz1DlLnW1eb4d7oyf3WVkIm5rpk6xbI60EfaoiFxWNfe+rGYfJKJjITmcZ7nob/qc/pPSH/4n3crXmo/WFP1N93/8S/NY3WD3l0gFZTNcWjHCgKb4Vos41DVbrpNHfYV/zRL5f+wSz1u/gx6d9lU0jpuCyG/z1xtO448nejRV1JxJp8zgnhMgG7geuAvKBrcCPpZQvpStPMaRUTs7e+JG6zxuv9qAedULXn0tBS7CFyubKdjGIaxXEO8xLRaG7kPE54xmfO54JuRPUOUedvQ4v2x5cxppVY5hY4ubHlfeyZP4SAL7z/77D/x3/AOVP+5k3YT8kEYlwMBKr+XcUgQAtlhBEwmaP3lMYgux8lxKBQhc5BW6yC93kFLmt6x64npg5k80/+A1bTvga59/evjL5vFtm8K/fb2TG1iXM/OmdsejR7hbZ5uu1QTd9PmSCQY9e90sNvQuBEFlZGN4sbFk9NOIx459w7/VieDwxo96XRFtF4x97DKBPusuOGWcWTJirDlD/u41Vlmh8APveUyJihqC+nNKD+6leU0DpvHq8m+6Byj+A3Ym3bi+lF0+j+nd/YdE170L9KHC4wW4dHa497dexZx6wu6xnLnUf/UwK55cep41JxV4mFXtTvl5rIJzQrdV+7XrLy9h59XgLA2BCIGxy3aPrWXHLHGYb22j9XB3+OiepU+8ZaWtJCCFeA04Gvgt8CtwILAAukVL+4yjS61VL4slv/R9GTSvg7Du+3fFByA9/v4s339rNQd80rp23G65eDtld+yvyhXwdjH+8GNT567rNT74rPyYEURGICkKOM6fLz7auW8/G7y/mg2kLmLGwiC/OmwfAa++uYdvyWsrKX2XsDdcSKh5Pc62f5vr21kBbc89bAU6PnRyr5p9dGN8S8JBTqNxSC2RsxokZCCKD0esAMsV97DoQ4GC9nTUVYzhh6xLGTi/GlpuL6fMRqtnPgQOSLdO/xkkHV1JQ+4ky8n5/j/PfJ9hsymB7POqclYXIsq49WbEwwwoTVrxQ9X5q//QnAEb98Id458yJEwAPwujfLsFjJdkYSl+Nq/Q7IT8c2ARV71O74kXcriq82dVJo7YedOKvc1JU1pL0+VFhOJKISaLouBKEJpU4JcQ7tB35xo9Yf/y93Lkul0MU8PIlJjPevQuuekx5gOiGjOxuEkJcCLwCXC6lXGmFCeAdoEhKWXYUafZKJN58+EF2bZzM1Fl72oWi5RA8tYA3PzbZ5buTqaOf5+wfPaL+IEBbuI19zfuStggOtx3u9jtznDmxFkBiiyDPlYeUkoi1OXwoGGnfDN46h4LtG8y3h0fYV1/N5s2rOK6qkMaCaRSO8dLWFKS1KQj0rA9fIPE4wmTZg3hsfrKEH49sxWM24w434w42YA+2YgbbDboMBDCDHa85xpkrFeO+SG5zBQUNu5I+r8+fSlPOBCbse73rhOz2doMdNegeD8KbYMw9HlUL93hUbb6DsY8TAqv7Rjh6P2Mq3phChtTAe8FADMAPOE37YcPj8PYD6v7EayB7JIQDav1TyA/huCPkV+HhQMdnoTbVSskAwnYvn3zhT2rMpRfjKpkqEn9GdTMVSinNuPBbgEeAE6SU23qZZq/HJDoIxZXz4clreG3vKPb67iRr7IuErzyDCqubaF9jFXUtDdhNB/aIA3fIiSvsxBVy4Aw7cUacuMIOsqSXQlsh+SKfbJGDBy8u6cEhXRAxCIclkTCEI+qIRARhUxAxBRHZfys8bWE/7kAdbn89rkAdbr+6dlvXzmAjhuxZt9KxIlwu63BiOF2xe8Pp7PDMbG7Gt/49AHIvvhjXlCmdDXeWt0NYVBCEMz0zzxIZ1DVwi8E4lbdbLDfnXPWYuo9e96Dm3QkzYolGQIlGB2Hx90B02lKIU6r0LHGSyRfjbTl3Gde/6e6xD6pMFYm1gJRSzksInwOsA66WUv4t4Vl31j8vLy+P3g5cv/nwg+z4uAxHcD8mDkKuEuzhBkAihRNpODENp3K9kU6kiS0SwGaGMCJBbGYQIxKyzu33YbuH2uKZAHxm74sU1W7FHajDHm5TbQohUhrldoPtVLNNnD18Fk3D6cJwxaXpcHa8d7l6VRPPuEHSXjIka+BDgXiBiIpCsrAMZ82uA3z3yfd48IppzDa3wHM3qQc3/J015vQeOyvMVJHYCeyUUl6cED4V2AncLqX8Q8KzfhEJQn6W3Lwcv+foHdoZkQC2SBDDDGGLGWx1tskwNiLqEBFswsRmmNiExG6T2GwSuw3sdrDbhTqiG8tbZ4fLhuF0KIMbXWjjdCAc1tlaeHOwzsbbb/qZ/sFvAdh22p2ce80ESo8vjIkCR9Fdkg4GexcNDNEa+GCnKzEYZEIRW0xnbEvaKlpjTh+8i+kskdghpbwkITwqErdJKf/YyzSPegrsE9/4Kc1CNWoKIh8zZlwRXk8WdqeBw2nH7rZhd9pxuB3Y3Q4cbjt2jxOHx4ndY9WWLYNtxK+UtB/7LmI9pWpHfafZP8lmCQ0GhkIXjSZDWf1rtX4ilQh8uko57xwE+3UDfdIqylSR6HV3Uw/SPCqRiB+XADoPZg8CEgUi3rAONqHQXTQaTQ/po1ZRdyKRro72rUCZEJ06+mda5y0DkYnEGU5n3/Ftps7aw66Nk3nz4QcHIgt9QtWqLZ0EAsB7+hxm/vROZmxdQtWqASnSY8a/ZXNKEfCePofSxYvxb9mchpxpNBlG9YbUIjDpP9Sz6g3H/DXpaklcBLwMfFlK+WJc+CpgpJSy125V+2QKbA+eZSK671uj0Rwtmbri+h/AW8ASIUQRajHdDcAZwJcGIgMHd9SnFIGz7/g2PPwgB3d0vxo6E+jO+HtPn6O7ZzQazVGRzhXXuSi3HFei3HJsQ7nleOEo00vPHtcajUYziMnUlgRSyibgDuvQaDQaTQaS2Q5jNBqNRpNWhpKrcBMQeXl56c6KRqPRDBoaGxtBLUlI2mgYSiIRRrWMut+EoTNRZWnsuxwNe3SZ9j26TPsWXZ6KXMCUUiYdfhgyInEsRF1+pBq40fQeXaZ9jy7TvkWXZ8/QYxIajUajSYkWCY1Go9GkRIuERqPRaFKiRUKj0Wg0KdEiodFoNJqUaJHQaDQaTUq0SGg0Go0mJXqdhEaj0WhSolsSGo1Go0mJFgmNRqPRpESLhEaj0WhSMqxFQgiRLYT4jRCiRgjRJoT4QAhxabrzNRgQQpwlhJApjuMT4p4rhFhnlfEhIcSfhBDD2l+OEGKsEOIhIcRqIUSLVW5npYh7nRDiYyGEXwhRJYR4QAjhThJvlBDicSHEESFEqxDiHSHEvH5/mQygp+UphChP8Zt9IEncYVue8QxrkQBWAguAHwAXoXbHWymEuDCtuRpc3APMTTjKow+tf9R/APuAS4DvAJcCrwghhvPvbwpwLdACvJEqkhBiIbAceBe4ALWb4zeBxxLiua10Pg98C7gMaAbeEELM6vvsZxw9Kk+LVXT+zf4uPoIuzziklMPyAC4EJHBZXJgAVgPb052/TD/nm9t3AAAGEklEQVSAs6zy+3I38d4DNgJGXNi51mevTvd7pLH84svjy1Z5nJUQxwbUAC8mhN9ixZ8TF3a7FXZyXJgL2Av8M93vmwnlaT0rB17oQXrDujzjj+Fck7sM5Uf+xWiAVL+Ex4HjhRDT05WxoYIQohQ4DVgqpTSj4VLK14Bq4Ip05S3dxJdHF5wOjEb9JuNZDoToWH6XAZullBviviMAPAmcK4TIObYcZzY9LM/eMKzLM57hLBIzgG1Jflyb4p5ruudPQoiwEKJRCPGyEOKUuGfRMtyS5HOb0WXcHUnLT0rpA/bQsfxmJMaz2IRqkZT1RwYHKWdb4xZBIcRmIcRtQgiREEeXp0XSnYiGCUXAziThdXHPNalpBH4NvI0qszLgXuBdIcTnpZTraS/DuiSfrwNOHoB8Dma6K7+ihLip4oH+PUd5GfgA1W1UBCwEfg8cB9wdF0+Xp8VwFglQfY5H82zYI6XciBpriPKOEOIlVO3rZ8AX46OnSqafsjfU6Gn56d9zN0gp70gIWimEWA7cKYT4tZSyIj56V0n1fe4yk+Hc3VRL8tpAoXVOVovQdIGU8gDwb1RfOqgyhtTlrMu4a3pTfvr3fPQ8jrKFs+PCdHlaDGeR2AqUJZmGOdM6J+uP1HSPQXsta6t1Tjb2MBNdxt2RtPyEEFnAZDqW39bEeBYzgQjwSX9kcIgQtQHx45O6PC2Gs0isBPJRc/fjuR7YIaXcNvBZGtwIIUajpreuA5BSVqH6fxfEi7EQ4hygFHg+HfkcRKwDDgCLEsKvBRx0LL+VwEwhxGejAUIIpxX3dSllUz/ndTBzPUog3o8L0+VpMZzHJP4BvAUsEUIUAZ8CNwBnAF9KZ8YGA1Y/7l5gA1APHI9aWOcBvhcX9R5UF9STQohHgDHAL4D1wDMDmedMQwhxpXV5mnX+vBCiGGiVUv5TShkWQtwLPCaEeBh4FjVB4BfAs1LKdXHJLUEtsnteCPE9VHfIXajy/soAvE7a6a48hRDXov63XwGqUF1HC1HrKn4lpayMS27Yl2eMdC/USOcB5AIPo2prfpTB63JxmD5iZXcv8BHQgJqzfwB4CpiRJO75KFHwA4eBPwMF6X6HdB+obrlkR3lCvIWoKcMB1PqSXwKeJOmNBpaiDJoPtTD0jHS/Z6aUJ2qs7HXUAsUgagX1GuCGFOkN6/KMHno/CY1Go9GkZDiPSWg0Go2mG7RIaDQajSYlWiQ0Go1GkxItEhqNRqNJiRYJjUaj0aREi4RGo9FoUqJFQqM5BoQQE63tL+9Ld140mv5gOK+41mg6IYTozcKhSf2WEY0mQ9AiodF0JNFP0pnA14FHgHcSnh1GrcT1AOH+z5pGM/BokdBo4pBSLou/F0LYUSKxNvFZHP5+z5hGkyb0mIRGcwwkG5OIDxNCfEUI8ZEQok0IsVsI8VUrznghxLNCiDohRLMQYlmyfZOFECVCiD8IISqt7Tb3CyEeEUKMHMDX1AxjdEtCo+k/LgZuRW2PWQd8DfiLECII3A+8Cfw3ymvpTagWyc3RDwshxgNrASfKK+keYApwG/AFIcSpUsrGAXsbzbBEi4RG03+UAdOltSWmEOJpYB/Ks+h3pJQPWvH+KIQoAK4XQvynlLLFCv8tat+IWVLtzYGVzjOovSbuBu4bkDfRDFt0d5NG03+8IOP2TJZSHgZ2oDa4+V1C3HdQgjARQAiRh2qJvAT4hRDF0QMoB3YD8/v7BTQa3ZLQaPqPvUnC6oEaKWUgSTi076s8DVWJ+5p19DR9jaZP0SKh0fQfkV6GA4iE8zLg8RRx244mUxpNb9AiodFkJrtRu6o5pZSvpzszmuGLHpPQaDIQKWUtah/2y4UQpyc+F4oRA58zzXBDtyQ0mszlNtS+yquEEE8AG1EVu88AXwKeQM9u0vQzWiQ0mgxFSrlPCHEKcA9KFBai1lLsA/4O/C2N2dMME4SUvfFnptFoNJrhhB6T0Gg0Gk1KtEhoNBqNJiVaJDQajUaTEi0SGo1Go0mJFgmNRqPRpESLhEaj0WhSokVCo9FoNCnRIqHRaDSalGiR0Gg0Gk1KtEhoNBqNJiX/H+35YX9EZe+kAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "catalysis_data.plot(x='Time', style='-x')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The theory of catalytic reactions guarantees that the total mass must be conserved.\n",
    "However, this is not the case in our dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    500.00\n",
       "1    415.09\n",
       "2    418.32\n",
       "3    494.71\n",
       "4    514.28\n",
       "5    565.40\n",
       "6    598.95\n",
       "dtype: float64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "catalysis_data.sum(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This inconsistency suggests the existence of an intermediate unobserved reaction product X.\n",
    "[(Katsounaros, 2012)](http://www.sciencedirect.com/science/article/pii/S0013468612005208) suggested that the following reaction path shown in the following figure.\n",
    "\n",
    "![](scheme.png \"Reaction Scheme\")\n",
    "\n",
    "The dynamical system associated with the reaction is:\n",
    "$$\n",
    "\\begin{array}{cc}\n",
    "\\frac{d \\left[\\mbox{NO}_3^-\\right]}{dt} &= -k_1\\left[\\mbox{NO}_3^-\\right], \\\\\n",
    "\\frac{d\\left[\\mbox{NO}_2^-\\right]}{dt} &= k_1\\left[\\mbox{NO}_3^-\\right] - (k_2 + k_4 +\n",
    "k_5)[\\mbox{NO}_2^-], \\\\\n",
    "\\frac{d \\left[\\mbox{X}\\right]}{dt} &= k_2 \\left[\\mbox{NO}_2^-\\right] - k_3 [X],\\\\\n",
    "\\frac{d \\left[\\mbox{N}_2\\right]}{dt} &= k_3 \\left[\\mbox{X}\\right], \\\\\n",
    "\\frac{d \\left[\\mbox{NH}_3\\right]}{dt} &= k_4 \\left[\\mbox{NO}_2^-\\right],\\\\\n",
    "\\frac{d \\left[\\mbox{N}_2O\\right]}{dt} &= k_5 \\left[\\mbox{NO}_2^-\\right],\n",
    "\\end{array}\n",
    "$$\n",
    "where $[\\cdot]$ denotes the concentration of a quantity, and\n",
    "$k_i > 0$, $i=1,...5$ are the *kinetic rate constants*.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "+ Assume that you are a chemical engineer and that you are assigned the task of designing a reactor for the conversion of nitrate to nitrogen. Before you start designing, you collect information in an attempt to characterize your state of knowledge about the problem. How many different sources of uncertainty can you think of?\n",
    "\n",
    "+ Which of these uncertainties would you characterize as aleatory?\n",
    "\n",
    "+ Which of these uncertainties would you characterize as as epistemic?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computational Model\n",
    "\n",
    "We will develop a generic computational model for the solution of dynamical systems and we will use it to study the catalysis problem. The code relies on the [Fourth-order Runge-Kutta method](https://en.wikipedia.org/wiki/Runge–Kutta_methods) and is a modified copy of [http://www.math-cs.gordon.edu/courses/ma342/python/diffeq.py](http://www.math-cs.gordon.edu/courses/ma342/python/diffeq.py) developed by Jonathan Senning. The code solves:\n",
    "\n",
    "$$\n",
    "\\begin{array}{ccc}\n",
    "\\dot{\\mathbf{y}} &=& f(\\mathbf{y}, t),\\\\\n",
    "\\mathbf{y}(0) &=& \\mathbf{y}_0.\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def rk45( f, y0, t, args=() ):\n",
    "    \"\"\"Fourth-order Runge-Kutta method with error estimate.\n",
    "\n",
    "    USAGE:\n",
    "        y = rk45(f, x0, t, args=())\n",
    "\n",
    "    INPUT:\n",
    "        f     - function of x and t equal to dx/dt.  x may be multivalued,\n",
    "                in which case it should a list or a NumPy array.  In this\n",
    "                case f must return a NumPy array with the same dimension\n",
    "                as x.\n",
    "        y0    - the initial condition(s).  Specifies the value of x when\n",
    "                t = t[0].  Can be either a scalar or a list or NumPy array\n",
    "                if a system of equations is being solved.\n",
    "        t     - list or NumPy array of t values to compute solution at.\n",
    "                t[0] is the the initial condition point, and the difference\n",
    "                h=t[i+1]-t[i] determines the step size h.\n",
    "        args  - any other parameters of the function f.\n",
    "\n",
    "    OUTPUT:\n",
    "        y     - NumPy array containing solution values corresponding to each\n",
    "                entry in t array.  If a system is being solved, x will be\n",
    "                an array of arrays.\n",
    "\n",
    "    NOTES:\n",
    "        This version is based on the algorithm presented in \"Numerical\n",
    "        Mathematics and Computing\" 6th Edition, by Cheney and Kincaid,\n",
    "        Brooks-Cole, 2008.\n",
    "    \"\"\"\n",
    "\n",
    "    # Coefficients used to compute the independent variable argument of f\n",
    "\n",
    "    c20  =   2.500000000000000e-01  #  1/4\n",
    "    c30  =   3.750000000000000e-01  #  3/8\n",
    "    c40  =   9.230769230769231e-01  #  12/13\n",
    "    c50  =   1.000000000000000e+00  #  1\n",
    "    c60  =   5.000000000000000e-01  #  1/2\n",
    "\n",
    "    # Coefficients used to compute the dependent variable argument of f\n",
    "\n",
    "    c21 =   2.500000000000000e-01  #  1/4\n",
    "    c31 =   9.375000000000000e-02  #  3/32\n",
    "    c32 =   2.812500000000000e-01  #  9/32\n",
    "    c41 =   8.793809740555303e-01  #  1932/2197\n",
    "    c42 =  -3.277196176604461e+00  # -7200/2197\n",
    "    c43 =   3.320892125625853e+00  #  7296/2197\n",
    "    c51 =   2.032407407407407e+00  #  439/216\n",
    "    c52 =  -8.000000000000000e+00  # -8\n",
    "    c53 =   7.173489278752436e+00  #  3680/513\n",
    "    c54 =  -2.058966861598441e-01  # -845/4104\n",
    "    c61 =  -2.962962962962963e-01  # -8/27\n",
    "    c62 =   2.000000000000000e+00  #  2\n",
    "    c63 =  -1.381676413255361e+00  # -3544/2565\n",
    "    c64 =   4.529727095516569e-01  #  1859/4104\n",
    "    c65 =  -2.750000000000000e-01  # -11/40\n",
    "\n",
    "    # Coefficients used to compute 4th order RK estimate\n",
    "\n",
    "    a1  =   1.157407407407407e-01  #  25/216\n",
    "    a2  =   0.000000000000000e-00  #  0\n",
    "    a3  =   5.489278752436647e-01  #  1408/2565\n",
    "    a4  =   5.353313840155945e-01  #  2197/4104\n",
    "    a5  =  -2.000000000000000e-01  # -1/5\n",
    "\n",
    "    b1  =   1.185185185185185e-01  #  16.0/135.0\n",
    "    b2  =   0.000000000000000e-00  #  0\n",
    "    b3  =   5.189863547758284e-01  #  6656.0/12825.0\n",
    "    b4  =   5.061314903420167e-01  #  28561.0/56430.0\n",
    "    b5  =  -1.800000000000000e-01  # -9.0/50.0\n",
    "    b6  =   3.636363636363636e-02  #  2.0/55.0\n",
    "\n",
    "    n = len( t )\n",
    "    y = np.array( [ y0 ] * n )\n",
    "    for i in range( n - 1 ):\n",
    "        h = t[i+1] - t[i]\n",
    "        k1 = h * f( y[i], t[i], *args )\n",
    "        k2 = h * f( y[i] + c21 * k1, t[i] + c20 * h, *args )\n",
    "        k3 = h * f( y[i] + c31 * k1 + c32 * k2, t[i] + c30 * h, *args )\n",
    "        k4 = h * f( y[i] + c41 * k1 + c42 * k2 + c43 * k3, t[i] + c40 * h, *args )\n",
    "        k5 = h * f( y[i] + c51 * k1 + c52 * k2 + c53 * k3 + c54 * k4, \\\n",
    "                        t[i] + h, *args )\n",
    "        k6 = h * f( \\\n",
    "            y[i] + c61 * k1 + c62 * k2 + c63 * k3 + c64 * k4 + c65 * k5, \\\n",
    "            t[i] + c60 * h, *args )\n",
    "\n",
    "        y[i+1] = y[i] + a1 * k1 + a3 * k3 + a4 * k4 + a5 * k5\n",
    "        y5 = y[i] + b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5 + b6 * k6\n",
    "\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Calibrating the Catalysis Model to the Experimental Data\n",
    "\n",
    "Now that we are certain that our generic ODE solver works, let us use it to develop a solver for the catalysis model. All, we need to do is define the right hand side of the dynamics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_catalysis(y, t, kappa):\n",
    "    rhs = np.zeros((6,))\n",
    "    rhs[0] = -kappa[0] * y[0]\n",
    "    rhs[1] = kappa[0] * y[0] - (kappa[1] + kappa[3] + kappa[4]) * y[1]\n",
    "    rhs[2] = kappa[1] * y[1] - kappa[2] * y[2]\n",
    "    rhs[3] = kappa[2] * y[2]\n",
    "    rhs[4] = kappa[3] * y[1]\n",
    "    rhs[5] = kappa[4] * y[1]\n",
    "    return rhs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try to calibrate the parameters of the model to the data, manually. Because the parameters are too small, let us work with the transformed version:\n",
    "\n",
    "$$\n",
    "\\xi_i = \\log\\left(180k_i\\right).\n",
    "$$\n",
    "\n",
    "Also, let's draw the graph corresponding to this model.\n",
    "We have the following variables:\n",
    "+ $\\xi$ corresponding to the scaled unknown parameters\n",
    "+ $y$ which is the prediction of our model at all timesteps for which we have data.\n",
    "+ $y_m$ which are the measured data.\n",
    "\n",
    "The graph will look as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: Catalysis Pages: 1 -->\n",
       "<svg width=\"134pt\" height=\"116pt\"\n",
       " viewBox=\"0.00 0.00 134.00 116.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 112)\">\n",
       "<title>Catalysis</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-112 130,-112 130,4 -4,4\"/>\n",
       "<!-- xi -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>xi</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"23.5\" y=\"-85.8\" font-family=\"Times,serif\" font-size=\"14.00\">ξ</text>\n",
       "</g>\n",
       "<!-- y -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>y</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"27\" cy=\"-18\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-13.8\" font-family=\"Times,serif\" font-size=\"14.00\">y</text>\n",
       "</g>\n",
       "<!-- xi&#45;&gt;y -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>xi&#45;&gt;y</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M27,-71.7C27,-63.98 27,-54.71 27,-46.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"30.5,-46.1 27,-36.1 23.5,-46.1 30.5,-46.1\"/>\n",
       "</g>\n",
       "<!-- ym -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>ym</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"99\" cy=\"-90\" rx=\"27\" ry=\"18\"/>\n",
       "<text text-anchor=\"start\" x=\"90.06\" y=\"-86.8\" font-family=\"Times,serif\" font-size=\"14.00\">y</text>\n",
       "<text text-anchor=\"start\" x=\"97.06\" y=\"-86.8\" font-family=\"Times,serif\" baseline-shift=\"sub\" font-size=\"14.00\">m</text>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.dot.Digraph at 0x1a26dd0990>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gc = Digraph('Catalysis')\n",
    "gc.node('xi', label='<&xi;>')\n",
    "gc.node('y')\n",
    "gc.node('ym', label='<y<sub>m</sub>>')\n",
    "gc.edge('xi', 'y')\n",
    "gc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "+ In the computational cell above add the edge that connects the model output the measurement. Make sure you understand why the arrow should point the way it does. Should it go from the model output to the measurement? Or should it go from the measurement to the model output?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6c2f795aca79485fb0e879be4b23852e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.359, description='xi1', max=2.0, min=-2.0, step=0.05), FloatSlider(v…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import interactive\n",
    "def compare_model_to_data(xi1 = 1.359, xi2 = 1.657, xi3 = 1.347, xi4 = -.162, xi5 = -1.009):\n",
    "    \"\"\"\n",
    "    Compare the model predictions to the data.\n",
    "    \"\"\"\n",
    "    t = np.linspace(0, 180, 100)\n",
    "    kappa = np.exp([xi1, xi2, xi3, xi4, xi5]) / 180.\n",
    "    y = rk45(f_catalysis, (500., 0., 0., 0., 0., 0.), t, args=(kappa,))\n",
    "    fig, ax = plt.subplots(figsize=(10, 10))\n",
    "    catalysis_data.plot(x='Time', ax=ax, style='s')\n",
    "    ax.plot(t, y[:, 0], color=sns.color_palette()[0], label='Model NO3-')\n",
    "    ax.plot(t, y[:, 1], color=sns.color_palette()[1], label='Model NO2-')\n",
    "    ax.plot(t, y[:, 2], color=sns.color_palette()[5], label='Model X')\n",
    "    ax.plot(t, y[:, 3], color=sns.color_palette()[2], label='Model N2')\n",
    "    ax.plot(t, y[:, 4], color=sns.color_palette()[3], label='Model NH3')\n",
    "    ax.plot(t, y[:, 5], color=sns.color_palette()[4], label='Model N2O')\n",
    "    plt.legend()\n",
    "    \n",
    "interactive(compare_model_to_data, xi1 = (-2, 2, 0.05), xi2 = (-2, 2, 0.05), xi3 = (-2, 2, 0.05),\n",
    "                                   xi4 = (-2, 2, 0.05), xi5 = (-2, 2, 0.05) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the calibration problem. \n",
    "\n",
    "### Questions\n",
    "\n",
    "+ Obviously, you do not want to be calibrating models by hand. Can you think of a \"natural\" way to calibrate a model?\n",
    "+ No matter what we do, we cannot really match the data to the model exactly. List at least two reasons why this is the case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uncertainty Propagation\n",
    "\n",
    "As discussed previously, there various reasons why a model cannot be calibrated perfectly. Some of these are:\n",
    "\n",
    "+ lack of data;\n",
    "+ the existence of measurement noise;\n",
    "+ the fact that the model is just not perfect.\n",
    "\n",
    "Ignoring for the moment the possibility that the model is just bluntly wrong, we see that the lack of data or the presence of noise will induce some uncertainty in the values of the calibrated parameters. We are going to represent uncertainty on parameters by assigning a probability density on them. There are systematic ways of estimating the uncertainty induced because of the calibration process, but this will not concern us now.\n",
    "For the moment, assume that somebody told us that the uncertainty in the scaled parameters $\\xi_i$ of the model is as follows:\n",
    "\n",
    "\n",
    "| Variable | Value            |\n",
    "|---------|------------------|\n",
    "| $\\xi_1$ |1.35 &plusmn; 0.05 |\n",
    "| $\\xi_2$ |1.65 &plusmn; 0.08   |\n",
    "| $\\xi_3$ |1.34 &plusmn; 0.11 |\n",
    "| $\\xi_4$ |-0.16 &plusmn; 0.16 |\n",
    "| $\\xi_5$ |-3.84 &plusmn; 0.20 |\n",
    "\n",
    "But what does this information actually mean? As we will discuss in the following lectures, this information can be used to assign a probability density on each one of these parameters, say $p(\\xi_i)$, that *models* our state of knowledge about them. For example, let us assume that our state of knowledge about $\\xi_1$ is given by a Gaussian probability density (don't worry about the notation - we will explain it in subsequent lectures):\n",
    "\n",
    "$$\n",
    "p(\\xi_1) = \\mathcal{N}(\\xi_1|\\mu_1=1.35, \\sigma^2 = 0.05^2),\n",
    "$$\n",
    "\n",
    "which we can visualize as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a262ed9d0>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEDCAYAAAAlRP8qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAdm0lEQVR4nO3de5BkZ33e8e9v7rNz3bt22WV1jyR0QSEKAgGWEbgKcZMQwVAoEAfb5ZSdUmKIiAEFQSzbJcDYchxsbKqEwWWKBElYBIwkyCbYICzFSEKsWAmtbstetDsz3TM7l52Z7jd/nHNmemZnpk/3nNPd77vPp2pqRqf79Hn7qPeZ37znfd9jzjlERMR/bc1ugIiIZEOBLiISCAW6iEggFOgiIoFQoIuIBKKjGQc1s3miXybjzTi+iIinBoGyc27F7LZmDFs0szJgQ0NDDT+2iIivisUigHPOrdi70pQKHRgfGhoaKhQKTTq8iIh/hoeHKRaLq/ZsqA9dRCQQqQLdzC43s3vM7JCZTZrZPjP7z2bWnXcDRUQknapdLmZ2AfB9YD/wH4DjwOuB24CLgPfl2UAREUknTR/6u4Ee4Abn3NPxtu+a2R7gPWb2AefcXG4tFBGRVNJ0uSRhXVy2vRg/Vsq0RSIiUpc0gf4lYBT4nJmdZWaDZvZ24P3AZ5xz5VxbKCIiqVTtcnHOPW9mVwL3AAcqHvo959wtK+1jZtXGI2oAukiDjE3O8sjBAq85dwud7RrYFrI0F0X3APcCR4DrgQLwC8DvmFl5tVAXkdbw4a89xn37jvLa87bwZze+gr7uZk0/kbyl+T/7B8AAcLlzbjrettfMAP6LmX3BOfds5Q7OueG1XjCu4FWli+RsrlTmvn1HAfjeU8f59S89zJc/8Erif78SmDR/f10O7KsI88TD8f4XZN4qEcnEvkNLJxX+w89GODZxskmtkbylCfRDwMVmtmHZ9lfF33+ebZNEJCsPPTsKwNaBbtriovzJoyea2CLJU5pAvwM4A/i2md1gZm8ws98FbgYecM79ONcWikjdkkB/zblb2LO5D4Anj040s0mSo6qB7py7B3gjcBL478DXiS6O/lfgulxbJyJ1c87x8LNjAFxx5ibO29YPwFMvKtBDlepyt3PuAeCBnNsiIhk6cHySkclZAK44cyOHCtPct++oulwCpkGpIoHafySqxPu62jlnaz/nbY8q9CePTtCM+yBI/hToIoE6VIgGpr1kYy9tbcb52wcAmJiZ5+i4RrqESIEuEqjDxRkAdgz1AnD21j7a46EuujAaJgW6SKAOF6MKfedwDwDdHe3s2RyNPlagh0mBLhKoQ4WlFTrAOVujfvTnRqaa0ibJlwJdJFBJhb5jqGdh2/bB6CZjL07MNKVNki8FukiA5kplXoyn+O8cXqzQtw1E4f6ipv8HSYEuEqAXJ06SjEysrNC3DcQVuka5BEmBLhKgw4XFtfQq+9C3xoF+7MRJjUUPkAJdJECH4iGLwxs66e1qX9iedLnMzpcZn55vStskPwp0kQAlFXpldQ6wLb4oCrowGiIFukiAkklFOyv6zwE293WR3NtCF0bDo0AXCdDCkMXhpYHe0d7G5j4NXQyVAl0kQEeKp04qSiQjXXTnovAo0EUCdPxEtGxuMqql0lYNXQyWAl0kQKPxOuib+7pOeWxhLLoq9OAo0EUCMz1bYnquBMDGlQJd0/+DpUAXCczo1OzCz5s2rFSha/p/qBToIoEZPVER6P2rd7noomh4FOgigUkq9M52Y6D71NsGJ10uEzPzTM+WGto2yZcCXSQwo5NR5b1xQxeWzCKqkIxDBxiZVJUeEgW6SGBGJ+cA2LTCBVGIgj5RmJprSJukMRToIoEZi4csrhboAz0dxLcWZaziAqr4T4EuEpiRONBXGrII0NZmDMdV+pgq9KAo0EUCs1ChrzBkMTG8oROAgir0oCjQRQIzWqXLBRb70ccmVaGHRIEuEphk2OLagR5X6NOq0EOiQBcJTJoKPelD1yiXsCjQRQJSKruFfvE0FbpGuYRFgS4SkPHpOcrxvZ/TVOga5RIWBbpIQJIhi5DuoqhGuYRFgS4SkMoulI1rDFtc6HKZVKCHRIEuEpCReKXFge4OujpW/+eddLmMz8wzXyo3pG2SPwW6SECK8TDEobgCX83GvsXHi9PqRw+FAl0kIEk4D/VWCfSK7hhdGA2HAl0kIGkDfbiigteF0XAo0EUCMj49D1QP9O6OdjZ0tQOq0EOSOtDN7Gozu8/MCmY2ZWb7zOzX82yciNQmbYUOFeu5qEIPRqpAN7P3Aw8ATwPvBt4K/Cmw+rgoEWm4JNAHUwS6VlwMz6k3HFzGzHYDnwM+4py7veKh7+TWKhGpS30VurpcQpGmQv9A/P1P8myIiKzf+Iwq9NNZmkB/HfAE8A4z229mJTM7aGZ/YGYrdrnE/eyrfgFDWb4JEYmMJ10uPVX/+F4IdI1DD0f1/+uwM/76E+AW4CfA64HfAXYD782tdSKSmnOupi6XwR4FemjSBHobMAC8xzn3lXjbXjPrBT5kZh93zv2scgfn3PBaL6gqXSR703Ml5krRUotpAj15TjLUUfyXpstlJP7+7WXbvxV//+fZNUdE6lUZzLUEuir0cKQJ9B+vst3i71rZR6QFVAZzmouigwr04KQJ9Lvi79cu234t4ICHMm2RiNSlMphrqdAnZuYoJ3fFEK9V7UN3zv2dmX0L+FMz28LiRdGbgD9zzj2XcxtFJIUk0Dd0tdPZXr1WSwK97ODE7PzCRVLxV5qLogD/CvgEcDOwFXge+Bhw+1o7iUjjjNcwwgVYEuDFqTkFegBSBbpzbhL4UPwlIi1oYdp/ymCuDP7i9By7c2mVNJJWWxQJRC1j0AEGejqweGhDMsNU/KZAFwlELQtzAbS1GQPd0R/p4xrpEgQFukggkio7bYUOGroYGgW6SCAW1nHpTTvWQZOLQqNAFwlErX3olc/V9P8wKNBFAlFPoGuBrrAo0EUCkVTZtYwnV5dLWBToIoGoq8slXhNdwxbDoEAXCcDsfJnpuRKwGNJpqEIPiwJdJACVFXZtfejRiBgFehgU6CIBWLJ0bg196IMLo1wU6CFQoIsEoNalc5c/d3x6Hue0hK7vFOgiAUgCvau9jZ7O9P+skwp9tlRmZk73qvGdAl0kAOMV67hYsuJWCstXXBS/KdBFAlDPtH9YGugauug/BbpIAOoZgw7LbnKhCt17CnSRANQb6F0dbfR2tkevMaVA950CXSQAybT/WgO9ch9V6P5ToIsEoNbbz1VaGLqoPnTvKdBFAlBvlwssXkhVhe4/BbpIANYT6OpyCYcCXSQA9dx+LqHb0IVDgS4SgGKd49Bhsd9ddy3ynwJdxHOlsmNiJr65xTq6XLRAl/8U6CKem6hz6dzl+6jLxX8KdBHPVXaVaNji6U2BLuK5JUvn1nC3ooQuioZDgS7iuSSI2wz6u2q/KJpU6FOzJeZKWkLXZwp0Ec8VK5bObWtLv3RuQkvohkOBLuK5pO+7nv5zWDrUUSNd/KZAF/HcemaJLt9PFbrfFOginltvoPd2ttPZbkteS/ykQBfx3HoD3cw0Fj0QCnQRz9V7+7lKC9P/ZzT932cKdBHPVY5yqdegpv8HQYEu4rnxdXa5VO6rLhe/1RXoZnarmTkzeyTrBolIbdbbhw6q0ENRc6Cb2cuADwNHs2+OiNQq6feudxw6wJDuWhSEmgLdzNqALwB/Cfw0lxaJSGrOuUwqdHW5hKHWCv0/AruAj+bQFhGp0eRsiVLZAQp0gdTjnMzsbOCTwHudc+Nmta8ZISLZqgzgdY1y6dESuiFIFegWpfdfAN92zt2T4vmFKk8ZSnNcEVlb5UXMTCr0KQW6z9JW6L8G/AvgohzbIiI1WlKh99Q/sSgJ9ImT85TLrq5VG6X5qn4CzGwLcDvw+8CkmQ1X7Nse//eMc24m2cc5N3zqKy15zQKq0kXWLQn0/u4OOtrrn1aSdNc4F4X6eqp9aZ40n4BdROH7+8BYxddVwMXxz7fm1D4RWcPCLNF1VOewtLtGY9H9leZT8DPgF1fY/kdAP/CrwPNZNkpE0hnPYNr/8v2L03PsXterSbNUDXTn3Alg7/LtyYVP59wpj4lIY2Qx7R9goLsDs6jLRUMX/aW1XEQ8lsWkIoC2NmOgO6rv1OXir7o73pxzV2fYDhGpQxYrLSaGNnQyPjOvCt1jqtBFPJas45LFqBTNFvWfAl3EY1l1ucDibFEFur8U6CIeyzLQk9fQ9H9/KdBFPFbM4PZzicUuF92GzlcKdBFPZbV0bmJQfejeU6CLeGpmrszsfBmA4Q1d6369Id21yHsKdBFPjU3NLvy8MYNA123o/KdAF/FUoWKp22ENWxQU6CLeKkxHFbpZNhOLkgW+itNzOOfW/XrSeAp0EU8lFfpgTyftGaxfnlTo82XH9Fxp3a8njadAF/FUEujDG7JZu3xo2YqL4h8FuoinkouiWYxwgVOX0BX/KNBFPJWEbhYXRGFphV7QvUW9pEAX8VRhoULPJtA729sYiC+Mjk3OVnm2tCIFuoinxuIqOosx6InktcZUoXtJgS7iqeJUdtP+Exv7kkBXhe4jBbqIp5LQ3ZhRl0vlaxUU6F5SoIt4qpBcFM2hy2V0Ul0uPlKgi3jIObfY5ZJphR4Fuip0PynQRTw0NVtithSttJjtRdHol4P60P2kQBfxUGE624W5Fl6rT6NcfKZAF/FQIeOlcxdfSxW6zxToIh5KZnK2GQuTgbKwKf7lUJyeo1TWiou+UaCLeKhQMQa9LYOVFhPJiBnntJ6LjxToIh7KemGuxMa+xf54dbv4R4Eu4qGkDz3LWaKwtD9eQxf9o0AX8dBIvHjWlv5sK/SeznZ6O9sBTS7ykQJdxEMjJ6JA39SXbaCDRrr4TIEu4qHRuELf3N+d+WsnC3Spy8U/CnQRDx0/cRKAzblU6Jpc5CsFuoiHFiv07AM9uWGGbnLhHwW6iGecc4uB3pd9l8smrYnuLQW6iGfGp+eZj2dx5nNRNFlCV4HuGwW6iGeOT55c+HlLDhdFtwxEr3n8hALdNwp0Ec+MVARt5czOrGyN++WPT5ys8kxpNQp0Ec+MxhX6QE8H3R3tmb9+UvVPnJxnZq6U+etLfhToIp5JukLy6G5Z/rrHVKV7pWqgm9k1Znanme03sykzO2hmd5nZJY1ooIgslVyszOOCKCz2ocPieHfxQ5oK/TeAlwKfBd4E/Hb83w+Z2ZU5tk1EVjCS46QigL6uxfVcdGHUL2lWxv9N59yLlRvM7D7gGeA/ATfk0TARWdlIjpOKAMyMLQNdvDA6rQrdM1Ur9OVhHm8rAE8Bu/JolIisLhnlksekokTSj66RLn6p66KomW0FLgYez7Y5IlJN3n3oUBHoqtC9UvPNCM3MgM8T/TL49CrPKVR5maFajysikZF42GJeXS5QGejqQ/dJPXeX/RRwHfArzrknMm6PiKyhVM53HZdEMrnomCp0r9QU6GZ2G/BB4Cbn3J2rPc85N1zldQqoShep2fETJ4mXcWH7YI596APqcvFR6j50M/sk8BHgZufcHfk1SURWc3R8ZuHn7UM9uR1HF0X9lCrQzezjwC3ALc65T+XbJBFZzZFiFOgbutoZ6K6nxzSdJNDHZ+Y5Oa/p/76o+okwsw8CtwLfAB5YNpnopHPuRzm1TUSWSSr07YM9ROMT8lF58+njJ2Z5yXBvbseS7KT5Ff/W+Ptb4q9KzwFnZtkgEVndkYVAz6//HJZN/584qUD3RNVAd85d3YB2iEgKR8ejPu0zBvPrPwcY6O6gt7Od6bnSkn57aW1abVHEI5VdLnkyM3YMR8c4XFSg+0KBLuKR5KJo3oEOsHMo6mY5VJzO/ViSDQW6iEeSCv2MHIcsJnbExzhcUIXuCwW6iCemZ0uMz8wD+V8UBdgRXwg9VFCF7gsFuognjlROKmpIl4v60H2jQBfxxJGKYN020IBAjyv0I+MzlJL1BqSlKdBFPPHiRBToW/q76OrI/5/uzniUS6nsdG9RTyjQRTyRVOiNqM4BdgwtTib6ufrRvaBAF/FEcnFyRwNGuAD0dXcw2BPNPTysoYteUKCLeOK50SkAXrp5Q8OOmfSja+iiHxToIp54fiQK9D2bGh/omlzkBwW6iAdKZccLY3Ggb+5r2HE1ucgvCnQRDxwuTjNXioYONqPL5WBhqmHHlPop0EU8kHS3mMGujY1byvbM+K+BZ45N4pzGorc6BbqIB5ILojuHeunuaG/Ycc/eGgX65GyJFzUWveUp0EU88Fxcob+0gRdEAc7a0kdyY6QDxyYbemypnQJdxAPPj0ZhuqeB/ecAPZ3tC8voHjh+oqHHltop0EU8kFTouxtcocNit4sq9NanQBdpcc65xTHoDa7QAc7ekgS6KvRWp0AXaXHHT8wycTJaB/3MBo5BT5y9tR+AA8dVobc6BbpIi9t3eByA9jbj3G39DT9+0uXywugUs/Plhh9f0lOgi7S4fYeiQD93az89nY0bsphIKvSyW7w4K61JgS7S4pIK/aKdg005/o7BHnrjXyRPHVU/eitToIu0uH2HigBctKM5gd7WZgu/TB49WGxKGyQdBbpIC5uanV+4GNmsCh3gsl3DADz6QqFpbZDqFOgiLWz/kQmSJVQubFKFDnDZ7iEAHjtY0P1FW5gCXaSFJf3nO4Z62NTX1bR2XL57IxCt6fK0xqO3LAW6SAtLujia1X+e2L2pd+EXyiPqdmlZCnSRFuWc4x9+NgLAvzxrU1PbYmZctivqdlGgty4FukiLeub4JD+Pbwz92vO2Nrk1cNnu6MLoPz031uSWyGoU6CIt6ntPHQdgS38XF5wx0OTWwKvO3gzAT49McHBMdzBqRQp0kRaVBPprzt1CW5s1uTXwij0bF/rR7993tMmtkZUo0EVa0FypzIMHov7z17RAdwtAR3sb11ywDVCgtyoFukgLun/fUU6cnMcMXnfelmY3Z8EvvewMAH74zCiFqdkmt0aWU6CLtKAvP/gcAK//Z9vYNtjT5NYseu15W+jtbKdUdtz72OFmN0eWUaCLtJinj53g+09H3S03Xrmnya1ZqqeznbdcugOAz//fp5kvaTndVqJAF2kxf/m9A0A0med157dG/3mlf3f1ObQZvDA6zd8+eqjZzZEKqQLdzPrN7A4zO2xm02b2sJm9Le/GiZxufnhghL/5xxcA+LdXnUV7C4xuWe7srf1ce0lUpf/RA08xPjPX5BZJIm2FfjfwXuBjwJuBfcDdZnZtXg0TOd2MTs7y4a89BsClu4b41y3W3VLppmvOo6u9jedHp7j5fzyGc1qwqxVUDfQ4tN8A/Kpz7gvOue8C7wd+AHwm5/aJnBZeGJ3iXX/+A54dmaKz3bj9nZfS0d66PaLnbR/gE29/GQB/95Mj3PSVR5iM73sqzdOR4jnXA0Xg68kG55wzsy8Cnzezi5xz+/JqoEioJmbm+PHBIt98/DBffeggs6UyXe1t3PGel3PBGc1djCuNd1+xm8d/XuSvf/g8f/voIX74zAg3vnIPv3jBNs7fPkBXR+v+QgqVVftTycx+QJThr162/ZXAg8AvO+e+uuyxaqv3DA0NDVEo1L7Izx/e/ySPHTx1v7XexmoP1ftn4mq7uVWPtMY+a7Z75Qfrea9rNK2+dq/+cque17X3qf04q+1Uz3Gi/eo433W0e65UZmxylpHJpeO4tw9289l3vZxXn9s6486rcc7xxe8/y23ffIK50tJ3PdTbyZb+Lvp7Omm36CbXyVebRV/WepcIGmKot5M/fvflNe83PDxMsVgsOueGV3o8TYW+GXhyhe2jFY83zI8PFti7/1gjDymSq0t3DfG2y3Zy45V7mnIT6PUwM/7NVWdxzYXb+asfPMs3HjvM4eIMAMXpOYrTumC6ki393bm8bppAhypFzykbVvntkYgr+KGUx17iTRfv4PztqyxUtMZve1vlwbUqhLWKh9X2W+04a+9T+4Ea1ba1jrXmPnWUXi3R7jo+J6u+1io7dbQZwxs62bWxl3O29jO8oXk3rsjK7k0b+OibL+Kjb76IQ4VpDo5Nc2ziJMcmZpiaK1EuO0plKJXLlFz08+l8IXVDV9rorU2aVx1h5So8WaB5dIXHcvOuK3Y38nAiUqOdw73sHO5tdjNOS2muWvwEuNDMlj/3kvj749k2SURE6pEm0O8GhoG3Ltv+PmC/RriIiLSGNF0u3wT+N/AFM9sMPEM0Dv01wNtzbJuIiNSgaqDHY86vA34v/hommin6DufcvTm3T0REUkp1qdU5Nw78VvwlIiItSFO5REQCUXWmaC4HNSsDNjRU11B0EZHTUrFYhKgnfMVivFmBPk/018F4HbsnvwWK2bUoeDpntdH5qo3OV23Wc74GgbJzbsXu8qYE+nok68RUm40qi3TOaqPzVRudr9rkeb7Uhy4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKB8G4cuoiIrEwVuohIIBToIiKBUKCLiASiZQLdzPrN7A4zO2xm02b2sJm9LeW+55jZPWZWNLMJM/ummV2Ud5ubrd5zZma3mplb4etII9rdDGa2y8z+2Mz+3sxOxO/36hr2f4WZfcfMJs1szMy+YmYvybHJTbWe82Vmd67y+Xow52Y3jZldE7/v/WY2ZWYHzewuM7uk+t7ZZVjLBDrRvUvfC3wMeDPRXZHuNrNr19rJzLYB3wPOJLo13nuATcD/MbNdeTa4BdR1ziq8EXhVxVfa/Xx0LtFn4wTwnVp2NLMLgb2AAe8Efg24HNhrZv3ZNrNl1H2+YidY+tl6FfCBzFrXen4DeCnwWeBNwG/H//2QmV251o6ZZphzrulfREHigOsrthnw98ATVfa9HZgGdlZs20y0NO/nmv3eWvSc3RrvO9zs99HA89VW8fN18fu/OuW+XwUOAX0V2y4ASsCHm/3eWvB83QkUmv0eGny+tq2wbRgYA75WZd/MMqxVKvTridYG/nqywUXv6ovABVX+9LgeuN85d6hi3xHgXuAd+TS3JaznnJ12nHPlevYzs07gLcD/dM5NVrzeT4EHgRuyaWFrqfd8na6ccy+usK0APAVUq7Izy7BWCfSLgX0rfIgeq3j8FGbWC5wDPL7Cw48B2+I/Z0JU1zlb5gkzK8V98H8R8Llaj7OBXlb/jKU5z6ejfjM7Gn++njOzzwTcPbUiM9tK9PlY6bOTPCfTDEt1k+gG2Aw8ucL20YrHV7KRqJthdIXHKvc95bdnAOo9ZwBPAx8BfgTMAlcBNwPXmNkrnHNjWTbUc8l5XO0z1mtmvc656Qa2qdU9CjxCFFLtRNdq/j3wWjO7yjk318zGNYKZGfB5oqL502s8NdMMa5VAh6iPrp7H1ruvz+p63865Ly3b9N14BMJ9wG8Cv5tB20Jzun7Gauac++yyTd82s/1EAffLwJcb36qG+xTRtYdfcc49keL5mXy+WqXLZYSVK8pN8feVfntBdMHB1bmv7+o9Zytyzt0PHCYajSCLRuLvq53raefcTAPb46svA2VOg8+Xmd0GfBC4yTl3Z5WnZ5phrRLoPwEuNLPl7UnGcK7YBxX/mXuAlfsxLwGOrXSxIhB1nbMq2oj+0cmiA0QjEFb7jNVznk9HFn8P+vNlZp8k6s682Tl3R7XnZ51hrRLodxMN8Xnrsu3vA/Y75/ZV2feNZnZGssHMNsWvdVfWDW0h6zlnpzCzXwK2E43ckFjc3/u/gBvMbEOy3czOJ6o2Q/6MZelGorwJ9vNlZh8HbgFucc59qoZdM8uwllhtMb6A8B3gUqKLc88QDbB/H/B259y98fP2Ar/gnLOKfbcTXYQ5BHwCmCeaaHM+cLlz7vnGvZPGWec5+xHwV8B+YA54NfAh4AhwRTzcKjhm9s74xyuIztmtRH/pTDrnvhU/51kA59yZFftdBPwjURh9GugDbgM6gZc75yYa8gYarJ7zZWZ7gC8Bf0N08b0deAPwW8D/A17nnJtv1HtoFDP7INFn4xtEn41KJ51zP4qft5c8M6zZA/IrBtIPAv+NKFRmgH8Crlv2nL3Ew62XbT+PaDz2ONEMtW8BL2v2e2rVc0b0j+0pYJJolMvTRDPcNjX7PeV8vtwqX89WPOfZyv+u2H4F8N34nBWIJhvtbvZ7arXzRTRq4654+3T8udwXB1Vvs99Tjudqb8rzlWuGtUSFLiIi69cqfegiIrJOCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQ/x/+KQqjvjh4pQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats\n",
    "from scipy.stats import norm\n",
    "xi1 = np.linspace(-0, 2, 200)\n",
    "plt.plot(xi1, norm.pdf(xi1, loc=1.35, scale=0.05))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that we do not believe that the value of the parameter can be less than 1.0 or greater than 1.6. Note that, we are deliberately trying to avoid the use of the term \"random\". There is nothing random in our example. Probability models a state of knowledge.\n",
    "\n",
    "How does this uncertainty propagate through the model? We will study this question with a simple numerical experiment. We are going to assign Gaussian probability densities on all the $\\xi_i$'s, sample them a few times, and run our catalysis model for each one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1c2e9e26159d403193f3f2910791ec21",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.359, description='mu1', max=2.0, min=-2.0, step=0.05), FloatSlider(v…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_samples(mu1 = 1.359, sig1=0.055,\n",
    "                mu2 = 1.657, sig2=0.086,\n",
    "                mu3 = 1.347, sig3=0.118,\n",
    "                mu4 = -.162, sig4=0.167,\n",
    "                mu5 = -1.009, sig5=0.368,\n",
    "                num_samples=1):\n",
    "    \"\"\"\n",
    "    Take a few samples of the model to study uncertainty propagation.\n",
    "    \"\"\"\n",
    "    fig, ax = plt.subplots(figsize=(10, 10))\n",
    "    catalysis_data.plot(x='Time', ax=ax, style='s')\n",
    "    t = np.linspace(0, 180, 100)\n",
    "    for i in range(num_samples):\n",
    "        xi1 = norm.rvs(loc=mu1, scale=sig1)\n",
    "        xi2 = norm.rvs(loc=mu2, scale=sig2)\n",
    "        xi3 = norm.rvs(loc=mu3, scale=sig3)\n",
    "        xi4 = norm.rvs(loc=mu4, scale=sig4)\n",
    "        xi5 = norm.rvs(loc=mu5, scale=sig5)\n",
    "        kappa = np.exp([xi1, xi2, xi3, xi4, xi5]) / 180.\n",
    "        y = rk45(f_catalysis, (500., 0., 0., 0., 0., 0.), t, args=(kappa,))\n",
    "        ax.plot(t, y[:, 0], linewidth=0.5, color=sns.color_palette()[0])\n",
    "        ax.plot(t, y[:, 1], linewidth=0.5, color=sns.color_palette()[1])\n",
    "        ax.plot(t, y[:, 2], linewidth=0.5, color=sns.color_palette()[5])\n",
    "        ax.plot(t, y[:, 3], linewidth=0.5, color=sns.color_palette()[2])\n",
    "        ax.plot(t, y[:, 4], linewidth=0.5, color=sns.color_palette()[3])\n",
    "        ax.plot(t, y[:, 5], linewidth=0.5, color=sns.color_palette()[4])\n",
    "    plt.legend()\n",
    "\n",
    "interactive(plot_samples, mu1 = (-2, 2, 0.05), sig1=(0.02, 0.4, 0.01),\n",
    "                                   mu2 = (-2, 2, 0.05), sig2=(0.02, 0.4, 0.01),\n",
    "                                   mu3 = (-2, 2, 0.05), sig3=(0.02, 0.4, 0.01),\n",
    "                                   mu4 = (-2, 2, 0.05), sig4=(0.02, 0.4, 0.01),\n",
    "                                   mu5 = (-2, 2, 0.05), sig5=(0.02, 0.4, 0.01),\n",
    "            num_samples=(1, 1100, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "Increase the number of samples from 1, to 10, to 100, to 1000. Each time you get a better description of uncertainty. This is a Monte Carlo simulation.\n",
    "\n",
    "+ Ok, the more samples you get the better your predictive error bars. But can you do this with any model? When would you face difficulties with such a program? What if you want to propagate uncertainties through a very complicated model, e.g., a climate model, which may take a few hours to complete a single simulation?\n",
    "+ Can you come up with any idea of accelerating the uncertainty propagation process?"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "number",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  },
  "widgets": {
   "state": {
    "0f9c50b41b864da1abfe0ff2fd704c74": {
     "views": [
      {
       "cell_index": 17
      }
     ]
    },
    "b95f55a04e394f1e9a18d215dfa8fc78": {
     "views": [
      {
       "cell_index": 22
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
